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EXECUTIVE SUMMARY 


A six degree of freedom (6DOF) computer simulation of the AIM-120 
AMRAAM has been developed to test the performance of various guidance laws using 
the kinematic boundary as a measure of effectiveness. Proportional navigation (PN) was 
used as the baseline for comparison. The effect of seeker noise on the PN law was 


studied. 


A velocity compensated PN law was tested against an angles only PN law and 
demonstrated that the velocity compensation will improve performance, but not to the 


level of the full PN law. 


A bang-bang law was tested as a continuation of earlier thesis work. This law 
performed poorly under the influence of drag, and would not be a candidate for use in a 


tactical missile. 


A modified PN law derived from differential games theory was tested that had 


lower performance than the PN law. 


An augmented PN law derived from optimal control theory was tested that had 
improved performance in the target's rear hemisphere and forward of 60 degrees relative 
to the nose of the target. This law did not improve the missile's performance against a 


cruise missile target. 


Preliminary work to extend the 6DOF simulation to include aerodynamic control 
of the missile was completed with the simulation capable of limited operation. More 


work needs to be accomplished to bring this model to full capability. 


XVII 


The 6DOF model was used to demonstrate the engagement of a theater ballistic 
missile by a RIM-67 STANDARD II (ER) missile. The STANDARD intercepted the 
target at a range of 2.2 meters off the nose, well within lethal range of the interceptor 


warhead. 
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I. INTRODUCTION 

The U.S. Navy’s experience with Japanese kamikaze attacks in the closing months 
of the Second World War demonstrated the woeful inadequacy of anti-aircraft artillery 
(AAA) against a massed modern air threat. Even with radar-controlled guns and the 
massed firepower of dozens of ships, the kamikaze were able to inflict heavy damage on 
the fleet. 

The advent of jet aircraft following the war exacerbated the aerial threat to surface 
units and changed forever the character of air-to-air combat. It was believed that the high 
speed and maneuverability of jet aircraft signalled the end of the dogfight and a 
requirement to engage targets at beyond visual ranges (BVR). The solution to both of 
these problems was some sort of guided missile. 

There are a number of ways to guide a missile so that it hits a target. The three 
simplest guidance laws are beam rider, pursuit, and proportional navigation. Beam rider 
guidance is most useful for a surface-to-air missile (SAM) installation. The launcher 
must keep the target locked in a radar beam throughout the engagement while the missile 
steers along, or rides, the beam. This requirement il] suits beam riders for the dynamic 
environment of aerial combat. Pursuit guidance requires the missile to turn so that it 
continuously points at the target. The missile may use some characteristic of the target 
such as its infrared (IR) signature, semi-active radar from the launching platform, or 
onboard radar to determine the targets relative position. As the name implies, this 
guidance law 1s most effective when attacking from the rear hemisphere of the target, or 
when attacking a stationary target, and it performs poorly in the targets forward 


hemisphere. 


Of the three basic guidance laws, proportional navigation (PN) is the most 
versatile, and therefore most frequently implemented. PN accelerates the missile laterally 
by an amount proportional to the angular rate of the line of sight from the missile to the 
target. PN or one of its various extensions or augmentations is the guidance law of 
choice in nearly all modern guided missiles. The reasons are simple. PN is: 

=> Cheap 

= Robust 

= Analytically tractable 

ТОШЕ ЕШ 

Optimal control theory promises to improve the performance of missile guidance 
systems. Optimal guidance laws, while the subject of extensive research, have yet to play 
a significant role in practical applications. Since optimal laws require an estimate of the 
target’s position, velocity, and acceleration (its state), they require computing horsepower 
that has only been available in miniaturized form since the late 1980’s. The computing 
requirements of optimal laws, and the successful extensions of PN laws have kept the 
optimal laws out of the mainstream, but modem, agile, stealthy aircraft and cruise 
missiles, and the growing need for theater ballistic missile defense (TBMD), have 
increased the interest in optimal guidance laws. 

The remainder of Chapter I examines the historical background, our goals in 
pursuing this line of research, and its benefits. Chapter II establishes the theoretical 
background for the simulation environments and for the various guidance laws we 
examined. Chapter HI descnbes our method of analysis using the kinematic boundary 


and our experimental procedures. Chapter IV presents experimental results and analysis. 


Chapter V presents our conclusions and suggestions for further research 1n this area. 


N 


А. HISTORICAL BACKGROUND 


Dunng the early 19505, the development of guided missiles was a major program 
for the U.S. military. SAM' developed during this era include the Army’s Nike family 
and Hawk, and the Navy's Terner, Tartar, Talos, and Standard. The primary air-to-air 
missiles (AAM) of the day were the Raytheon Sparrow, developed for the Navy, and the 
Hughes Aircraft Falcon, developed for the Air Force. Both of these systems were 
complex radar-guided missiles (Falcon had an IR variant) and suffered from many 


developmental problems that would be familiar to systems engineers today. 


While the engineers at Raytheon and Hughes were overcoming their technical 
challenges, a small team of scientists and engineers at the Naval Ordnance Test Station 
(NOTS) in China Lake, California, now the Naval Air Warfare Center Weapons Division 
(NAWCWPNS) began work on what would become one of the most successful AAM’s in 
history. Sidewinder (AIM-9) began as an after work project on a non-existent and 
frequently purloined budget with no official standing [1]. In the view of the air power 
theorists of the day, the age of the dogfight was over, so why would there be a need for a 
short range dogfight missile? The Vietnam War would soon prove the theorists wrong 


and demonstrate the value of Sidewinder. 


Sidewinder was designed from the beginning to be simple, reliable, rugged, and, 
above all, inexpensive. The motor, warhead, and fins were adapted from a stock five 
inch High Performance Air Ground (HPAG) rocket. The fins were modified with a 
mechanical device called a "rolleron" which minimized the missile’s roll rate without the 


need for additional electronics [1]. Most of the design effort went into the guidance and 


control section which was bolted on to the HPAG rocket as a unit, and Incorporated 
several innovations, including: 
e Torque balance servo control fins which provided the commanded control 
forces at all altitudes without complex electronics 
e Single gyroscope seeker which integrated the IR sensor and directional gyro 


e IR aiming reticles which reduced the missile’s tendency to guide on the sun or 
clouds 


The design was so simple that NOTS technicians would tell Air Force and 
Hughes personnel that the only test equipment they required was a flashlight and a 
Simpson meter [1]. While this may have been an exaggeration for psychological effect, it 
was not far from the truth. The first production Sidewinders cost the government $2,400 
and by the third year of production, the price was down to $1,400 per missile [1]. 
Тодау 5 advanced Sidewinders cost in the tens of thousands of dollars. Compare this to 
over $300,000 for an AIM-120 AMRAAM. Sidewinder was such a successful design 
that it was copied wholesale by the Soviet Union as the K-13 (NATO AA-2 Atoll), and 


used as the basis for the Israeli Python [2], [1]. 


Sidewinder' guidance law is a form of PN using only line of sight angular rate 
and a fixed navigation constant or gain. This guidance law is suitable for a dogfight 
missile with a range on the order of 5.5 km (18,000 ft.), but not for longer ranges. The 
general PN law incorporates the missiles closing velocity with the target in the 
computation of the gain and must be provided a measurement of the range rate. This is 


the realm of the radar-guided missile. 


The first radar-guided missiles in the U.S. inventory were the Air Force’s Falcon 


(AIM-4) and the Navy’s Sparrow (AIM-7). Both missiles used semi-active radar homing 


(SARH) seekers. The launch aircraft must illuminate the target throughout the 
engagement for these missiles to guide successfully. Doppler processing of the 
illuminator’s return from the target aboard the missile provides an estimate of the closing 
velocity. The need to continuously illuminate the target means that the launch aircraft 
must continue to close with the target during the engagement. This creates an obvious 


problem if the target s weapons have similar ranges to those of the launch aircraft. 


Falcon enjoyed a long career, retiring in 1988. Sparrow is still in use today, and 
as NATO Sea Sparrow is the point defense missile system aboard many U.S. and NATO 


ships. 


During the 1960's, Hughes began development of the missile that eventually 
became the AIM-54 Phoenix. Phoenix includes a strapdown inertial measurement unit 
(IMU) that allows its autopilot to steer the missile on course with periodic updating from 
a SARH seeker. In the terminal phase, the missile switches to an onboard active pulse 
Doppler radar. Finally, the missile has a simple data link with the AWG-9 radar aboard 
the F-14 launch aircraft that allows the aircrew to command the missile to perform 
several functions. All of these improvements permit the F-14 to simultaneously guide six 


missiles to different targets up to 176 km (110 miles) away. 


Raytheon’s AIM-120 Advanced Medium Air-to-Air Missile (AMRAAM) is the 
current generation of missile technology in the U.S. inventory. AMRAAM incorporates 
an IMU, a data link, and a pulse Doppler terminal seeker. Because its data link is more 
sophisticated than Phoenix, there is no need for a SARH seeker. In certain scenarios, 


AMRAAM is truly a "fire and forget" missile, using its IMU to fly to a point where the 


active seeker can take over. Generally, AMRAAM is launched with an initial intercept 
solution programmed into the autopilot by the aircraft radar and mission computer. Once 
fired, the data link can update the autopilot with target position while the launch aircraft 
turns away or engages other targets. Once the terminal guidance seeker is activated, the 
missile 1s completely autonomous. AMRAAM has substantial onboard computer 
processing available and can employ advanced signal processing algorithms and guidance 


laws. 


Proportional navigation can be shown to be an optimal solution under a set of 
limited conditions. Chief among these limitations is the assumption that the target does 
not maneuver during the engagement. This is clearly unrealistic, and there have been 
many extensions to the basic PN law to counter this limitation. Optimal control theory 
makes it possible to account for target maneuvers in the guidance law. This requires an | 
estimate of at least the target's acceleration and in some cases the complete target state. 
A range of tracking filters including the Alpha-Beta-Gamma and Kalman filters is 
available to provide these estimates. Single chip microprocessors and digital signal 
processors have made it possible to implement these guidance laws in the limited volume 
of a missile’s guidance section. Despite these developments and the potential advantages 
of optimal guidance laws, the practitioners have been slow to implement new designs. 
Some of this lag is due to the successful extension of the PN law, but much is due to the 
aversion of more experienced engineers for abandoning a technique that works in favor 


of techniques that have yet to prove themselves [3]. 


Modern agile aircraft like the MiG-29 and stealthy aircraft like the F-117 and F- 


22 may in some cases be able to defeat AAM’s using PN laws. It is thought that optimal 
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and hybrid guidance laws may overcome the limitations of PN laws. Optimal laws may 
also increase the range at which cruise missiles can be engaged, and developments in 
differential games theory (a field of mathematical optimization) may help solve the TBM 
problem [4]. 


B. GOALS AND BENEFITS 


The research presented in this thesis was motivated by three primary goals. First 
to create a set of 6DOF models for evaluating missile guidance laws, second to explore 
the use of the kinematic boundary as a measure of effectiveness (MOE) for evaluating the 
performance of the simulated missiles, particularly to compare optimal guidance laws 
with PN laws, and third to demonstrate an application of the models to a TBM 


interceptor. 


Much of the literature in the missile guidance field involves the use of two- 
dimensional simulations. While such models are fairly simple to set up and analyze, and 
are nol as computationally intensive as 6DOF models, they have difficulty simulating the 
effects of drag and aerodynamic contro! forces on the missile. Our goal was to create a 
simplified 6DOF model for guidance law development and testing, and a full 
aerodynamic model that would simulate both the aerodynamic control forces and the drag 
forces acung on the missile. The modular design of the Simulink® models makes it 
possible to test not only guidance laws, but autopilots, thrust profiles, and the effect of 


noise anywhere in the system on performance. 


There are a number of ways to construct MOE's for the evaluation of a missile's 
performance. Controls engineers would compute a cost function based on the miss 


distance, control effort, and possibly time of intercept. While the number produced by 
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such a cost function is useful as a basis of comparison, to the layman it is simply a 
number. For the warfighter, the engagement envelope is of paramount importance. The 
kinematic boundary represents the maximum range at which the missile will achieve a hit 
when there is no noise in the system. It is a graphical representation of which guidance 
law has the best performance. If several points in the boundary are tested using noise, the 
mean effect of the noise can be calculated and its effect on the engagement envelope 
demonstrated. This information can then be used to determine if one guidance law is 


truly more effective than another. We have used the kinematic boundary as the MOE 
throughout the AAM simulations. 

The final goal of this research was to provide a missile simulator that could be 
used in other research conducted for Navy TENCAP (Tactical Exploitation of National 


Capabilities) in the TBMD field. 


I. BACKGROUND 


A. SIX DEGREE OF FREEDOM (6DOF) DYNAMICS 


Newtons laws for both translation and rotation describe the motion of a body in 
three-dimensional space. There are three axes for translation, x, y, and z, and three axes 
for rotation, longitudinal, lateral, and vertical, giving nse to displacement in roll, pitch, 
and yaw respectively. These are the six degrees of freedom. The coordinate frame for 
these dynamics is centered on the aircraft center of gravity (c.g.) and fixed to the airframe 
with the x-axis on the nose, y-axis on the right wing, and z-axis pointing down. It is 
called the aircraft-body centered or ABC frame. This is a rotating frame in inertial space 
and for objects in different ABC frames to interact; their motion must be transformed into 


an inertial frame. 


For short ranges (« 200 km) the North-East-Down, or NED, frame is suitable. 
This frame assumes a flat earth, and reasonable altitudes so that gravity is a constant. A 
NED has its x-axis pointing north, y-axis pointing east, and its z-axis pointing down 
toward the center of the earth. An aircraft headed north in level flight will have pitch roll 
and yaw angles of zero degrees. A z-axis which points down seems counter-intuitive at 
first, but makes sense when one considers that this allows right hand turns to have ап 
increasing heading as seen on a compass. We will use the NED or flat earth 


approximation for the air-to-air engagement simulations. 


If the NED coordinate system were placed on the surface of the earth, it would 
become a rotating frame with the earth’s angular velocity. For long ranges and ballistic 
missile work, one final translation to the earth-centered inertial or ECI frame is required. 


In this fixed frame the x-axis points at the vernal equinox or first point in Aries (which is 
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really in Pisces), the y-axis is 90 degrees to the east, and the z-axis extends through the 
North Pole. We will use the ECI frame for the TBM interceptor demonstration. Figure 
2.1 shows the relationships of the three coordinate frames. 

2 


ECI Coordinate Frame 


Note: (0,0,0) in NED is 


ye the tangency point of the 
a nom NED plane and the Earth 
ABC Coordinate Frame / пету еј 
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Геше sls Relationship of ABC, NED, and ECI Coordinate Frames. 


There are four vector equations which describe the dynamics of a body in three- 
dimensional space [5]. They are the force equation, the moment equation, the attitude 
equation, and the navigation equation. The equations shown below are for the flat earth 


approximation. The individual terms are defined in List of Symbols and Abbreviations. 
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F 
Y, z—0,v,tB.2, Tw ( force) 


ој = ж. T pa (moment) (2.1) 
4 а 75 5224 (attitude) 
Рмер = B; v B (navigation) 


The attitude equation can be computed using quaternions as shown here or using 
Euler angles. The Euler angle formulation involves a singularity in the rotation matrix 
(Вв) when the missile passes through the vertical that does not occur in the quaternion 
formulation. Since the STANDARD missile is fired from a vertical attitude, the 


quaternion formulation will be used throughout. 


For the TBM interceptor demonstration, the round earth equations shown below 
in state space form are used. Note the addition of terms using 42r, which is the cross 
product matrix accounting for the Earth's rotation and the B matrix instead of Bg that 
rotates the ABC frame to NED coordinates, and then to ECI coordinates. Definitions of 


the individual terms are listed in the List of Symbols and Abbreviations. 


T 
P CN ( В r) 0 0 p 0 
у, : = BO% = О, + BQ -B 7 0 у, у Bg(p) | 2) 
à, 0 0 —J e : 0, Jen i 
д 0 0 02 -29|4 0 


These equations assume constant mass and a fixed center of gravity. Simulation 
of a missile that burns fuel and has a shifting c.g. as a result involves the addition of 
terms to the force and moment equations. For simplicity we have assumed a constant 


mass missile. 
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The MatLab@ functions FLATEARTHDYN.M and SIXDOFDYN.M implement 


these equations in the SimuLink® models described below. 


The equations of motion assume the motion takes place ina vacuum. As a result, 
there is no direct coupling between the force equation and the moment equation. A stable 
missile body with its c.g. forward of its center of pressure (c.p.) tends to act like a 
weather vane and align itself with the relative wind. In the simplified 6DOF model this is 
modeled by feeding back the angle of attack, which is the angle between the missile body 
and the velocity vector, and its derivative as a moment that steers the missile into the 
relative wind. The specifics of this feedback will be outlined below. The full 
aerodynamic model does not require this feedback as it generates the normal forces on 
the missile by generating a moment using control deflections and using the subsequent 
change in angle of attack to generate the forces. 


В. MISSILE MODELING 


The simulation environments are capable of modeling any missile the researcher 
chooses to represent. For this research, the AIM-120 AMRAAM, and RIM-67(ER) 
STANDARD II (SM-2) missiles were chosen. These weapons represent today's front 


line U.S. Navy technology. 


The model dimensions have been simplified to comply with the supersonic 
aerodynamic models in Zarchan, and Blakelock, but are generally representative of the 
actual missiles [7], [8]. The performance specifications are also simplified and based on 
capabilities reported in the open source literature, and on engineering approximations. 


They are in no way intended to be representative of the actual capabilities of these 


missiles. No official use U.S. Government or contractor proprietary documentation was 
used in the establishment of the model performance parameters. 

1. Airframes 

AMRAAM is a conventional missile design with fixed stub wings mounted 
forward on the missile body and controllable tail fins mounted aft. There are four wings 
and four fins mounted at 90-degree intervals around the missile body. Figure 2.2 shows 
the overall plan view of the missile, and the MatLab® file MISSILEDAT A.M establishes 
the model’s dimensions as required by Zarchan. The definitions of the dimensions used 


in Zarchan’s equations are shown in Figure 2.3 [7]. 


Missile Plan View AIM-120 AMRAAM 
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Missile Plan View RIM-67 STANDARD Missile 
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Figure 2.2. | AMRAAM and STANDARD models. Drawings to scale for comparison. 
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Figure 2.3. Dimensions and forces on a tail-controlled missile. From [7]. 


For the computation of the moments of inertia, the missiles are modeled as thin 
rods for the y and z-axes, and cylinders about the x-axis. The thin rod model was chosen, 
because the fins are not major contributors to the moment of inertia about the axes 
normal to the longitudinal axis, and the missile is much longer than its diameter so the 
contribution of the radius for the cylindrical model is minimal. The cylindrical model 
was chosen for the longitudinal axis because there is no moment of inertia for an 
infinitely thin rod about the longitudinal axis. Since the missiles are symmetrical, there 
are no cross terms in the inertial matrices. 

The model for SM-2 is more complicated. The extended range version of the 


missile is equipped with a large booster with controllable tail fins. Figure 2.4 shows the 
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SM-2(ER) model with the booster attached. Note that the winss and tail fins for the 
missile forward of the booster have been modeled as a single wing with a length equal to 
the wing plus tail fin and an area equal to wing plus tail fin. MISSILEDATA4.M 


contains the dimensions for the missile in this configuration. 


SM-2(ER) with booster attached 
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Effective Center of Pressure 





meters 


Figure 2.4. SM-2(ER) model with booster attached. 


Once the booster stage falls away, the SM-2 looks like the second drawing in 
Figure 2.2. MISSILEDATA3.M contains the dimensions for the missile in this 
configuration. For comparison, line drawings of the actual missiles are shown in Figure 


29 
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AIM-120 AMRAAM 


RIM-67 STANDARD II 
А 





Figure 2.5. AMRAAM and SM-2 (ER). 


2 Propulsion 

Both AMRAAM and STANDARD use solid fuel rocket motors. The actual 
missiles use dual propellant grain motors that provide a relatively high value of thrust 
initially to accelerate the missile to speed quickly, and then a lower level of thrust to 
sustain speed throughout flight. For simplicity, the motors are modeled as single grain 
motors of intermediate thrust values. 

Solid fuel motors used in military missiles must have a Department of Defense 
(DoD) Hazard Classification of 1.1 or 1.3 for use aboard ship [9]. According to Sutton, 
typical fuels of this type have specific impulses in a range of 180-270 seconds [9]. The 
thrust F produced by a rocket motor is given by:[9] 

F =I mg. (2.3) 
This equation assumes a constant propellant mass flow rate throughout the motor run. 

Assuming propellant mass fractions of 50 percent for AMRAAM and SM-2 
without its booster, and 80 percent for the SM-2 booster, with a six second burn time for 
AMRAAM and 10 seconds each for SM-2 and its booster yields the data shown in Table 


2.1. The thrust values chosen for use in the simulations are within the range of 
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feasibility, and were chosen to accelerate the missiles to their maximum speed in а 


reasonable time. 


(180<[5<270) Thrust (N) 
AMRAAM 23,062 — 34,594 23,000 


STANDARD II 62,209 — 93,314 80,000 
SM-2 BOOSTER | 137,655 — 206,482 180,000 


Table 2.1. Missile Thrust Values. 







3. Aerodynamics 

a. Simplified 6DOF Model 

The aerodynamics for the simplified 6DOF simulation are modeled as a 
feedback path from the missile state vector. The ABC velocities are used to compute the 
pitch and yaw angles of attack (a and D) that are then differentiated and fed back as a 
proportional-differential (PD) controller to the torque input of the missile dynamics block 
(See Appendix B, Thesisl.mdl) This feedback loop models the missile's natural 
tendency to act like a weather vane when the lift and side forces change the velocity 
vector and hence the relative wind. The lift and side forces are generated by multiplying 
the guidance law command accelerations by the missile's mass. 

The angle of attack response of the missile to a step input is similar to a 
second order response with a damped oscillation. This is also similar to the response of 
the full aerodynamic model to a step input on the control fins. The feedback gains were 
chosen to give the missiles a settling time of approximately 2.5 seconds, or an 


approximate first order time constant of 0.5 seconds. 
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Drag is modeled with two components, parasitic drag, that due to the 
misssile's shape and cross section, and induced drag, that caused by the generation of lift 
and side (normal) forces. The drag force D along the velocity vector is computed using 


the following equation [10]. 


2 


V 
D - (C4, mo E Ó REF (2.4) 


~ 


Since the steady state angles of attack generated by this model are small, 
less than one degree, the small angle approximation has been used and the cosine of the 
angle of attack has been set to one for computing the component of drag along the x-axis 
of the missile. 

Cao 1s computed using typical values provided in [6]. The data were faired 
to a polynomial curve using MatLab®, and the function DRAGTHESIS.M is used to 
compute the parasitic drag in the model. Figure 2.6 shows the variation of Cao with Mach - 
number. The upper curve is the result of the increased drag caused by turbulence around 


the missile's tail when the thrust plume is absent. 
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Parasitic Drag Coefficient (Ca) vs. Mach number 
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Figure 2.6. — Variation of parasitic drag coefficient with Mach number. 


Cy, 1s computed in two regimes, subsonic, and supersonic. Normally, Cg 
is a function of angle of attack, but in this simplified model, the angle of attack values are 
not realistic, and therefore, a different approach is required. 

For subsonic flight, Cg; is computed as the applied normal force in 25 
times the maximum value of Cgo in subsonic flight. This crude approximation only 
affects the missile for very short periods of time as it is subsonic only at launch and 


perhaps at the very end of an engagement. 
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In supersonic flight, a more accurate approximation based on the normal 
forces is used. The normal force coefficient Cy is computed using the following 
equation: 


Fy 


— (2.5) 
p V "A REP 


С, =2 
where Fy is the applied normal force. Cy; is then computed using the following [10]. 


Cy 
Cr = ze(lAR) (2.6) 


Since there are normal forces on both the y and z-axes, Cai is computed for 
each axis and the results are added to produce the value of Ca; used in Equation 2.4 
above. Figure 2.7 shows the parasitic and induced drag forces on the AMRAAM model 


for a missile in level flight executing a variety of turns at load factors up to 30 25. 
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Figure 2.7. Drag forces on the AMRAAM mode! for various load factors. 


b. Full Aerodynamic Model 

The full aerodynamic model follows the development in Zarchan for 
generation of both the aerodynamic moments and forces. Moments are generated by the 
deflection of the appropriate control surfaces (rudder or elevator). The simulated missile 
flies in a vertical attitude with the elevator surface horizontal and the rudder surface 
vertical. AMRAAM flies in a "cross" configuration with the tail surfaces at 45 degree 
angles to the vertical for ease of loading and carriage aboard aircraft. Modeling this 


involves a more complicated autopilot and the change to a vertical attitude does not 


materially affect the simulation. The aerodynamic moment T caused by a control surface 
деПеспоп 15 отуеп ђу:[7] 


om о Зик (2.7) 


Cy is a function of the angle of attack and the control deflection and is 


given by:[7] 
ТУ 0 
См = 204 cx — X cos) (X co -X crs ) 
= REF (2.8) 
So S+ lG + Ó 
n (Х со – Херу )+8 - (Ха-Хш) 


S REF p S REF 





Where p is a normalized speed for supersonic travel given by:[7] 


B - мае 1 (2.9) 


The normal force Fy on a body is given by:[7] 
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V 
РЕ = Cy PS rer (2.10) 


Cy; is again a function of the angle of attack and the control deflection and 


is given bv:[7] 
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The equations for Cy and Cy given above are valid for the supersonic 
regime. No such approximation based on missile dimensions exists for the subsonic 
regime. For subsonic speeds, the coefficients normally are determined empirically using 
wind tunnel or computed fluid dynamics data. Since these data were not available, Cy 


апа Су 1п the subsonic range are modeled as linear functions of the angle of attack and 


tO 
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control deflection [10]. Equations 2.7 and 2.10 are then used to generate the moments 
and forces. The values chosen for the coefficients are therefore arbitrary, but as with the 
drag model above, since the simulation spends very little time in the subsonic regime, the 
effects of this approximation will be minimal. 

The drag model is quite different from the simplified 6DOF model. 
Parasitic drag is computed in the same fashion as above. Subsonic induced drag is 
computed using Equation 2.6, because the model now explicitly calculates Cy. 


Supersonic induced drag follows an approximation given in Anderson [10]. 


C F] (2.12) 


Since the angles of attack are generally greater than one degree, the 
induced drag force due to each normal force is computed separately, and its component 
along the longitudinal axis of the missile is computed before being added to the other 
component. The parasitic drag force is multiplied by the cosines of both angles of attack 
to determine its longitudinal component. 

According to Stevens and Lewis, once the moments and forces have been 
determined, the 6DOF equations are solved in the following order [5]: 


e Force and moment equations 
Attitude equation 
e Navigation equation 


4. Guidance, Navigation, and Control 
a. Guidance 
Guidance laws are implemented as Matlab@ functions which compute the 
command accelerations (ne) for both lateral and vertical guidance. The inputs to the 
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guidance law are provided by the seeker head, which computes target range, range rate, 
azimuth, elevation, and angular rates from the actual target and missile state vectors. 
Measurement noise can then be added to any of the six output channels to study its effect 
on guidance law performance 

Guidance laws for the simplified 6DOF model must also generate the 
applied force on each axis for the computation of the drag forces. Simulink® generates 
“algebraic loop" errors when the forces are fed back from the input of the "Missile 
Dynamics" block (Figure B.1). Guidance laws for the full aerodynamic model do not 
require this additional output. 

Guidance laws requiring a tracking filter incorporate the filter's estimate of 
the target state and missile body frame accelerations as additional inputs. 

b. Navigation 

The inertial measuring unit (IMU), air data computer (ADC), 
accelerometers, and rate gyros provide navigation data to the missile simulation in the 
form of Euler angles, missile total velocity, acceleration, position, angles of attack, and 
body axis rotation rates. The IMU is mounted at the missile's c.g., thus simplifying the 
calculation of the accelerometer data. Although this research assumed a noise-free 
navigation system, noise sources could be added to any of the output channels to study 
the effect on performance. In particular, the effect of navigation system noise on the 
tracking filter could be studied. 

C. Control 

The simplified 6DOF model does not require an autopilot, since the 
guidance law command accelerations are directly converted into aerodynamic forces. For 
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the full aerodynamic model, it is necessary to convert the command accelerations into 
control deflection angles. For this purpose, an autopilot for tail-controlled missiles 
presented in Blakelock was adapted for use [8]. Blakelock’s autopilot contains feedback 
loops for a missile which does not guide during boost, and to correct for accelerometers 


which are not at the c.g. These loops have been deleted in this model. 


C SIMULATION ENVIRONMENTS 


Two distinct simulation environments were developed for this research. The 
simplified 6 DOF model was designed initially for the purpose of developing and testing 
guidance laws prior to using them in the full aerodynamic model. Problems with the non- 
linearity of the full aerodynamic model delayed its completion, and as a result, most of 
the simulation results presented were obtained from the simplified 6DOF models. All 
simulations operate in continuous time using the Simulink® ode45 Dormand-Price 


differential equation solver. 


The 6DOF models, THESISI.MDL (Figure A.1) and THESISIFILT (Figure 
A.14) employ the flat earth approximation (Equation 2.1) for their missile dynamics, and 
are streamlined models providing only the minimum number of subsystems required to 


quickly test guidance law operation. 


All the air-to-air simulations use a point mass target simulation developed in [6]. 


The target dynamics are modeled with the following vector equation:[6] 
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The target s lateral accelerations are modeled as a turn rate, œ, while the vertical 


acceleration is an input to the subsystem, a;. 


The TBMD model, THESISTBM.MDL (Figure A.23) uses the spherical earth 
model (Equation 2.2) for its missile dynamics. The target model used in this simulation 
involves a six dimensional state vector to simulate the dynamics of a point mass ballistic 


missile with no drag as shown below. 
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The full aerodynamic model, THESIS3.MDL (Figure A.18) presented the greatest 
design challenge. In order to meet Stevens’ and Lewis’ requirement that the 6DOF 
equations be solved in the proper order, the flat earth dynamics block was completely 
redesigned (Figure A.21) The moments and forces on the missile are computed as 
outlined above, and then fed to the missile dynamics block as inputs. This should have 


resulted in a model that could be run in both open loop and closed loop operations. 


Unfortunately, 1t was not possible to successfully close the loop with either the autopilot 


adapted from Blakelock, or any of several other autopilot designs. 


It was possible to control the missile laterally, and for short periods vertically in 
an open loop by using the control deflection angles as inputs. It is likely that the failure 
of the closed loop operations was due to the inherent non-linearity of the model, and 
possibly the order in which Simulink® solves the computations in the various Matlab® 
function blocks. It may be possible to code both the aerodynamics and missile dynamics 
blocks as one inline Matlab® function to overcome this failure. This model is presented 
here as a point of departure for future research. 


D. GUIDANCE LAWS 


Five guidance laws were examined during this research, proportional navigation 
(PN), velocity compensated proportional navigation (VCPN), bang-bang, differential 
games (DG), and augmented proportional navigation (APN). The PN laws were used to 
establish baseline performance for comparison with the other guidance laws. The bang- 
bang and VCPN laws were examined as an extension of thesis work by Swee [11]. The 
DG and APN laws are derived in the optimal control literature and are the focus of using 


the kinematic boundary as measure of effectiveness [12], [13]. 


The geometry of a typical air-to-air missile engagement is shown in Figure 2.8. 
The object of the exercise is to steer the missile using only lateral accelerations in such a 
way that it hits the target. The steering commands should be optimal in some sense, 


minimizing miss distance at least, and possibly control effort (divert) or time to intercept. 
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Figure 2.8. ‘Typical missile engagement geometry. From [13]. 


ДЕ Proportional Navigation 

PN provides steering commands to the missile, which are proportional to the 
angular rate of the target’s line of sight relative to a fixed reference. The command 
acceleration n, 1s given by:[7] 


a: 


coso, 


(2.139 


n. 


The cosine term in the denominator corrects the acceleration from the line of sight to the 
missile’s y-axis. 
PN with N=3 has been shown to be optimal and guarantees a hit under the 
following conditions:[14] 
e non-maneuvering target (no drag) 


• missile speed greater than target speed 
e target remains in missiles forward hemisphere 
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If the value of N' is sufficiently large, PN will always intercept a maneuvering 
target under these conditions. Blakelock implements PN as a turn rate, but his 
recommended values for a navigation constant equate to values of N' between three and 
five [8]. Higher values produce little improvement in performance. One other 
shortcoming of PN is that it does not account for the effect of the missiles dynamics 
(time constant) on the navigation solution. 

2: Velocity Compensated Proportional Navigation 

VCPN is an attempt to extend the basic PN law and account for the effect of drag 
on the missile. By adding a compensation term related to the missile's deceleration and 
the line of sight angle, the effect of the drag on the line of sight rate can be reduced. The 
VCPN law is given by:[13] 


А 


COSO, 


n -V, tano, (2.16) 


In his thesis, Swee showed that if the range rate information VC is available to the 
missile, VCPN is no better than the basic PN law [11]. 

S. Bang-bang 

Bang-bang guidance is a modification of PN in which the missile applies its full 
acceleration in the direction of the rate of change of the line of sight. The controls 
essentially "bang" on their stops whenever they are applied. This law would be useful in 
missiles that are controlled by thrusters that are not throttled and are either on or off. The 
bang-bang law is given by:[12] 


onlV.o 
п = А sgn(Vc 6) (2.17) 
COS O, 


The bang-bang law use here is modified slightly because of the effects of drag. 
First, there is a dead band of 0.01 degrees per second in the line of sight rate before the 
guidance law takes effect, and second the acceleration at ranges greater than 5 km from 
the target is restricted to 5 g’s. Inside 5 km, the acceleration is 30 g’s. This was done to 
prevent the missile from expending all of its thrust overcoming the drag from 30 g turns 
immediately after launch. 

4. Differential Games 

Bryson and Ho develop a guidance law based on differential games theory in 
which the pursuer (missile) seeks to minimize a cost function based on the miss distance 
and the control effort while the evader (target) seeks to maximize the cost function and 
thus survive. Both players are assumed to have perfect knowledge of the other's state. 
Under these conditions, the evader S optimal strategy is to match the pursuer turn for turn 


as shown here [12]. 


estos) (2.18) 





The pursuer’s optimal strategy is more complicated, but after assuming the 
pursuer can turn at a faster rate than the evader, that minimum miss distance is infinitely 
more important than minimum control effort, and linearizing about a nominal collision 


course, the resulting control law 15:[12] 


n ms V. б (2.19) 





This is a variation of PN where the navigation constant can be varied statically at 


launch based on an estimate of the target’s ability to maneuver, or dynamically with a 
30 


real-time estimate of the targets acceleration. Bryson and Ho say that c, and с, аге 
constants related to the respective energies of the evader and pursuer, but closer analysis 
shows that they are also related to the ability to maneuver or available acceleration [12]. 
The control law implemented here uses a value of 30 g for c, and estimates c, from the 
output of the tracking filter. 

5. Augmented Proportional Navigation 

This guidance law is drawn from Lin, Reference [13], and is a simplification of an 
optimal guidance law that accounts for both target maneuver and missile dynamics. The 
APN law used here does not account for missile dynamics. It uses a twelve-dimensional 
state vector with the targets relative position, relative inertial velocity, target inertial 
accelerations, and missile body frame accelerations. The tracking filter estimates the first 
three, and the body frame accelerations are provided by the accelerometers. The 


guidance law is given by the following vector equation:[13] 





The navigation constant A is computed for the full optimal guidance law as a 
function of f,,,. and the weighting functions on the miss distance and control effort. When 
the product of the weighting functions approaches zero, the navigation constant is equal 
to three. We have chosen a value of five to be consistent with the baseline PN law. 

The position and velocity components of the state vector are relative to the missile 


and in inertial coordinates; therefore, they can be computed directly from the seeker 


3] 


ranges and bearings. The time to go, to, is computed from the seeker range and range 
rate estimates. 

Only the y- and z-components of the control are used. The x-component is 
ignored. Note that in this form, the missile accelerations are not used. A constant 
diagonal matrix is used in place of the "0" matrix to add the effect of the missile time 


constants in the full optimal guidance law. 


E. TRACKING FILTER 


The tracking filter is based on an alpha-beta-gamma filter design by Bar-Shalom 
and Li [15]. This is a constant gain filter and therefore it is less computationally 
intensive than an adaptive filter like the Kalman filter. Zarchan recommends the use of 
constant gain filters, in part because of computational load and also because of stability 
[7]. The DG and APN guidance laws require this tracking filter for their estimates of the 


target's acceleration 


The filter is implemented as a MatLab® function ABGFILTER.M. It is 
interesting to note the use of the global variable XLAST to preserve the state estimate 
from time step to time step. The values for the filter gain were chose by trial and error 
from nomograms in Reference [15] to give the filter an initial settling time of less than 
two seconds as these simulations are initially noise free, and the choice of gains is 
dependent on the characteristics of the noise. The filter is a discrete time filter with a 
sampling frequency of 10 Hz. This was accomplished by placing the filter block between 
two zero order hold blocks. A sample of the filter's estimate of target acceleration with a 
6 g turn three seconds prior to intercept is shown in Figure 2.9. 
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Figure 2.9. o-fB-y Filter Performance. 
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III. GUIDANCE LAW TESTING 


А. KINEMATIC BOUNDARY 


Kinematics is the branch of mechanics dealing with pure motion without 
reference to the masses or forces involved. For the purposes of this research, a missile’s 
kinematic boundary is the locus of points representing the maximum range at which a 
target may be successfully engaged as a function of relative bearing from the target at the 
start of the engagement given a noise-free guidance and control system. To the pilot, this 
is the "firing envelope,” a critically important piece of information as it determines not 
only the success of an engagement, but the tactics required to prosecute the target. We 
chose the kinematic boundary as our measure of effectiveness for this reason. To the 
warfighter, graphs of average miss distance or control effort may be meaningful if he 1s a 
controls engineer, but a comparison of two guidance laws showing one to have a 
significantly larger firing envelope is far more useful. Figure 3.1 below shows a generic 
kinematic boundary (a circle) and is representative of the plots used in Chapter IV. The 
azimuth angles represent the relative bearing of the shooter from the target at the start of 


the engagement. 


А successful engagement has a miss distance of less than 5 meters for these 
simulations. This figure is based on the warhead of the AMRAAM having a lethal radius 
of approximately 10 meters, and the size of a typical modern jet aircraft. Figure 3.2 
shows the relationship of the 5 meter radius to a MiG-29 fighter. Clearly, a warhead 
exploding within 5 meters of the MiG-29 will do substantial if not fatal damage to the 


ай сга, 
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The Matlab® program files KBOUTER2.M and KBFILTER.M generate the 
kinematic boundaries. The resolution in range is 10 meters, and in azimuth is 5 degrees. 
These values were chosen as a compromise between speed of execution and plot detail. 
At these resolutions, a kinematic boundary can be generated in 9-12 hours with a 


Pentium® III, 700 MHz processor. One-degree resolution requires 48-60 hours, and 1 


meter would take approximately 8-10 times longer. 


A Representative Kinematic Boundary 
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Kinematic boundary. The shooter is on the boundary pointing at the target 


Figure 3.1. 
at the start of the engagement. 
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5 meter lethal radius 
around c.g. of MiG-29 





Figure 3.2. | Comparison of a 5 meter warhead lethal radius to a MiG-29 aircraft. 
MiG-29 drawing is from [2]. 


B. TEST SCENARIOS 


Candidate guidance laws are tested in three engagement scenarios: 
e Non-maneuvering target, co-altitude at 6,000 meters 


e Non-maneuvering cruise missile target at 50 meters 
e Maneuvering target, co-altitude at 6,000 meters 


By 


The non-maneuvering target engagement is used as a baseline for comparison of 
performance. The engagement begins with target and shooter at 6,000 meters altitude, 
approximately 20,000 feet, and Mach 0.83. These values would be typical of an intruder 


making a high altitude ingress to a target, and a combat air patrol (CAP) on station. 


The cruise missile engagement is intended to demonstrate the interceptor’s ability 
to engage a low-altitude non-maneuvering target like the Tomahawk missile. The 
AMRAAM is launched from the CAP station at the target, which is at 50 meters, 


approximately 150 feet. 


The maneuvering target engagement is the true test of missile performance. In 
this scenario, the target initiates a 6 g turn or "jink" toward the missile three seconds prior 
to impact. This turn toward the missile is most advantageous to the target as it forces the 
missile to make a tighter turn and expend more energy to keep up with the target than a 
turn away. The timing was chosen for two reasons. First, given that the missile will 
activate its terminal radar between 5-7 seconds pnor to impact, and the time required for 
the targets sensors to detect the radar, alert the pilot, and have the pilot take evasive 
action, the aircraft would be established in its maneuver about three seconds prior to 
impact. Secondly, the missiles settling time 1s modeled to be 2.5 seconds, so a maneuver 


at three seconds puts increased stress on the guidance law to keep up with the maneuver. 
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С. CANDIDATE GUIDANCE LAWS 


Table 3.1 shows how the guidance laws were tested. 


co-altitude cruise missile co-altitude 
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Table 3.1. Guidance Law Test Plan 
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The VCPN and bang-bang laws were tested to confirm earlier work by Swee in 
his thesis [11]. 


D. NOISESTUDY 


A study of the effect of seeker noise on missile performance was conducted using 
the PN (N =5) guidance law. The study was run at the 135-degree azimuth test point with 


100 realizations. The standard deviations of the noise signals were as follows: 


e Range 50 meters 

e Closing velocity 2 meter/second 

e Bearing I degree 

e Bearing rate .01 degree/second 


39 


THIS PAGE INTENTIONALLY LEFT BLANK 


40 


IV. COMPARISON AND ANALYSIS 


A. PROPORTIONAL NAVIGATION LAWS 


The PN guidance laws were tested to provide a baseline for comparison with the 
other guidance laws. Figure 4.1 shows the kinematic boundaries for the three PN laws 
against a non-maneuvering co-altitude target. Figure 4.2 shows the kinematic boundaries 
against the co-altitude, maneuvering target described above. Figures 4.3 and 4.4 are 


amplifications of the differences in performance of the three laws. 


Тһе N=3 guidance law is the poorest performer of the three. While N23 has 
been shown to be optimal for a non-maneuvering target, the effect of drag on the missile 
is similar to a target maneuver along the line of sight. The discontinuities or "divots" in 
the N —3 boundary are caused by drag slowing the missile more rapidly on those attack 
azimuths than others resulting in the missile slowing below the targets speed and 


stopping the simulation. 


Clearly, the N=5 law is an improvement for both scenarios. There is a slight 
improvement between N 25 and N=7 with a mean value of 315 meters for the non- 
maneuvering case, and 1,749 meters for the maneuvering case. The improvement from 
N=3 to N=5 has a mean value of 2,076 meters, non-maneuvering, and 7,555 meters 
maneuvering. Because of the relatively poor performance of the N 23 law, N=5 will be 


used as the comparison baseline for the other guidance laws. 


Figure 4.5 is a comparison of the performance of the N’=5 law against a co- 
altitude target and against a cruise missile target at 50 meters altitude. The omni- 
directional reduction in range is due mainly to increased drag as the missile descends into 


the heavier air at lower altitudes. 
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Figure 4.1. Kinematic boundary comparison of proportional navigation laws vs. non- 


maneuvering, co-altitude target at 6,000 meters and Mach 0.83. 
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Figure 4.2. Kinematic boundary comparison of proportional navigation laws vs. 


maneuvering, co-altitude target at 6,000 meters and Mach 0.83. Target maneuver was a 6 
g turn toward the missile at /,,73 seconds. 
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Comparison of PN Performance vs. a Non-maneuvering Target 
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Figure 4.3. PN law performance vs. non-maneuvering target. 
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Figure 4.4. PN law performance vs. maneuvering target. 
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Kinematic Boundary Comparison 
PN (N'25) vs Non-maneuvering and Cruise Missile Targets 
90 
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—— Non-maneuvering 


180 


Figure 4.5. | Kinematic boundary comparison of PN vs. non-maneuvering and cruise 
missile targets. 
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В. VELOCITY COMPENSATED PN LAWS 


The VCPN laws were tested as a continuation of Swee's thesis research [11]. A 
PN law with a fixed navigation constant, angles only, no Vc information, was tested 
against the non-maneuvering, co-altitude target. This simulates a guidance law like that 
used by Sidewinder. The gain was computed with N=5 and a fixed closing velocity of 
750 meters per second. VCPN laws with and without Vc information are compared to 


this law, and to PN with N=5. 


Figure 4.6 shows the kinematic boundaries for each of these laws. The VCPN 
law without Vc information is clearly an improvement over the angles only PN law, while 
the addition of Vc information to the VCPN law actually reduces the range. Since the 
incorporation of Vc information also includes the deceleration of the missile along the 
line of sight, the velocity compensation term adds nothing to the guidance laws 


performance. Neither of the VCPN laws performed as well as the full PN law. 
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Kinematic boundary comparison of VCPN laws vs. a non-maneuvering, 


Figure 4.6. 


co-altitude target at 6,000 meters and Mach 0.83. 
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e BANG-BANG 


The bang-bang guidance law was also tested as a continuation of Swee’s 
work.[11] The bang-bang law is used throughout the engagement to determine the effect 
of drag on its performance. As seen in Figure 47, bang-bang is clearly outperformed by 
the baseline PN law. Because of the aerodynamic drag on the missile, the guidance law 
must expend more energy in the end game when the line of sight angular rates begin to 
increase. There is a synergistic effect: as the angular rate increases, the missile must turn 


harder, generating more drag, which causes the angular rate to increase. 


Note that the bang-bang laws performance is highly aspect dependent. The 
enhancement in the target s forward hemisphere is most noticeable. The effect could be 
useful in TBMD work where the goal is to place the interceptor ahead of the target. 
Further, for an exo-atmospheric interception, the effect of drag on the bang-bang law 


would be greatly reduced. 
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Figure 4.7. 
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р. DIFFERENTIAL GAMES 


The differential games law was tested against both the non-maneuvering and 
maneuvering co-altitude targets. Figure 4.8 shows its performance against the non- 
maneuvering target, and Figure 4.9 against the maneuvering target. In both cases the 


performance showed no improvement over the baseline PN law. 


This law is a modification of PN with scheduling of the navigation constant based 
on the tracking filter’s estimate of the target’s total acceleration. It is clear that gain 


scheduling is not sufficient to increase the kinematic boundary of the PN law. 
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Figure 4.8. | Kinematic boundary comparison of the differential games law vs. a non- 


maneuvering, co-altitude target at 6,000 meters and Mach 0.83. 
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Figure 4.9. | Kinematic boundary comparison of the differential games law vs. 


maneuvering, co-altitude target at 6,000 meters and Mach 0.83. Target maneuver was a 6 
g turn toward the missile at t,,=3 seconds. 


Е. AUGMENTED PROPORTIONAL NAVIGATION 


The APN law was tested against both the non-maneuvering and maneuvering, co- 
altitude targets. Against the non-maneuvering target, the APN law’s performance is 
identical to the baseline PN law except for a small "divot" at 100 degrees. Figures 4.10 
shows the results against the maneuvering target. The jagged boundary is an artifact of 
the azimuthal resolution, and is smoothed out when the resolution is reduced to one 


degree. 


The APN law 1s clearly better in the targets rear hemisphere and forward of 60 
degrees relative to the nose. In the forward quarter from 90 degrees to 60 degrees there is 
a reduction in performance compared to the PN law. The mean improvement in the APN 
law 1s 4.45 km for all aspects, 1.24 km in the forward 120 degrees, and 8.48 km in the 
rear hemisphere. 

Figure 4.11 shows the results of APN against the cruise missile target. The 
kinematic boundary for the APN law is clearly smaller than PN law. Azimuth resolution 
was reduced to 10 degrees for this comparison to keep the APN simulation under 48 


hours in real time. 


9 | 


5) 


' APN (A= 
РМ (№ 


Н 
і 
5 
i 
1 


5) 


150000 


5) and PN (N' 
90 
270 


— 
c 


Kinematic Boundary Comparison 
АРМ (А 





180 


Figure 4.10. Kinematic boundary comparison of APN vs. maneuvering, co-altitude 
target at 6,000 meters and Mach 0.83. Target maneuver was a 6 g turn toward the missile 
at 1,573 seconds. 
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Figure 4.11. 


Kinematic Boundary Comparison 
APN and PN (N’=5) 
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Kinematic boundary of APN and PN vs. cruise missile target. Azimuth 
resolution is 10 degrees. 


F. NOISE STUDY 


A study of the effects of noise on missile performance was conducted at the 135- 
degree azimuth test point. Using the range for the kinematic boundary of the baseline PN 
law, it was not possible to hit the target in 100 realizations. The test point was moved in 
approximately 2 km to 45,500 meters and another 100 realizations were generated. 
Figure 4.12 is a scatter plot of the x and y miss distances for the 100 realizations. Figure 
4.13 1s the distribution of the Euclidean miss distances. 92 percent of the samples were 
within the required miss distance of 5 meters to be called hits. 
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Figure 4.12. Scatter plot of x and y miss distances for a noisy seeker. 
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Noise Study Miss Distance Histogram 
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Figure 4.13. Histogram of missile miss distances with a noisy seeker. 100 realizations. 
Probability of hit is 92%. 
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С. TBMD DEMONSTRATION 
Figure 4.14 is a demonstration of the 6DOF models ability to simulate a TBM 


interceptor. The target missile was launched from the equator on a northeasterly heading. 
The range of this missile is approximately 400 km. The interceptor was launched from a 
position 150 km north of the target launch site. The target’s initial velocity vector, [vy vy 
v;] in ECI coordinates, was [1200 10 1000]. The velocity profiles for the target and 
interceptor are shown in Figure 4.15. The interceptor was steered toward the target's 
apogee for the first 30 seconds of flight, and then followed the baseline PN law to an 
interception 2.2 meters ahead of the target. The plot is in ECI coordinates, the surface of 


the earth is approximately the bottom grid. North is to the nght. 
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Figure 4.14.  TBMD engagement by RIM-67 SM-2 (ER). Miss distance at intercept 
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Figure 4.15. Interceptor and target velocity profiles for TBM demonstration. 
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V. CONCLUSIONS AND FUTURE RESEARCH 


A. CONCLUSIONS 


The kinematic boundary is a natural, intuitive method of comparing the 
performance of guidance laws. Its form is immediately recognizable to the warfighter, 
and provides exactly the information required to make an informed decision as to which 


guidance laws would be of operational value. 


The VCPN laws showed the expected improvement over the fixed gain PN law, 
and the VCPN law with range rate information did not perform as well as the VCPN law 
without range rate information. Neither law performed as well as the baseline PN law. 
The bang-bang law showed an unusual kinematic boundary with ranges in the target’s 
forward hemisphere greatly extended over the rear hemisphere. The aspect dependence 
of this law and approximately 50 percent reduction in range throughout most of the 


envelope make this law a poor choice for guidance from launch to intercept. 


Under the conditions simulated here, an optimal control law, the augmented 
proportional navigation law, will perform better than a proportional navigation law 
throughout most of the kinematic boundary. Overall, the APN law's kinematic boundary 
was 4.45 km better than the PN law, on average. In the forward 120 degrees, the average 
improvement was 1.24 km and in the rear hemisphere, it was 8.48 km. This represents a 
1.4 percent improvement over the PN law for head-on engagements, and a 25 percent 


improvement for rear hemisphere engagements. 


The 6DOF model has been demonstrated as a test platform for evaluating 


guidance laws for use in the TBMD arena. 
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В. FUTURE RESEARCH 


This work suggests several lines of future research. First, the guidance laws 
tested here should be tested with a noisy seeker to determine the effect of noise on their 
performance. Second, the full aerodynamic model designed here needs to be taken to the 
point where it is fully operational. Third, the TBM simulation could be used to study the 
effects of guidance law selection on the Navy’s Linebacker TBMD capability. Fourth, 
the models could be used for a comparison of the kinematic boundaries of other missiles 
systems, particularly those that are potentially hostile to look for possible tactical 
advantages. Finally, the models could be used to test new guidance laws that will be 


developed future researchers. 


APPENDIX A. SIMULINK® MODELS 


The block diagrams in this appendix represent the four models used in this 


research. Sub-blocks which are not changed from the simplified 6DOF model in later 


models are not included with those models. 


pages: 


eos Simpliiicd ODOR... л 
• 6DOF with tracking filter......... 
e Full aerodynamic model........... 


e TBMD interceptor model ......... 


The four models begin on the following 
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Simplified 6DOF Aerodynamic Force Generator. 
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Figure A.2. 
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Simplified 6DOF Drag Model. Function blocks are DRAGINDUCED.M 


Figure A.3. 


and DRAGTHESIS.M. 
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Flat Earth 6DOF missile dynamics. 
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Figure A.4. 
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Internal missile dynamic model. SIXDOFDYN.M is the function block. 


Figure A.5. 
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Aerodynamic moment feedback. 


Figure A.6. 
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Missile IMU and air data computer. Function blocks are Q2EULER.M 
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and ALPHABETA.M. 


Figure A.7. 
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Missile seeker model. 


Figure A.8. 
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Figure A.9. 
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Figure A.10. Range, range rate, and time-to-go. Function block is TGO.M. 
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Target dynamics model. Function block is DYNAMIC3D.M. 


Figure A.11. 


74 


€G:91 000Zz-d3S-Z| pəluud 


| L/} гбед 





| рыл }8\вәц | \ўзом\гсдееш\.ә 


} 5 ——— — —_—_ 





о Y: РАИ 


rr — "nt 0 € p e аА а а 


0J9Z 0] Jas 'ulni ои 20у 
риоэәѕ јад зитреј ш ӘДЕ ШІП) 
јабје) оу) 5195 »20ја 5141 


ayey uinj 
YOUMS 








eyes usin} 


pjoysasyy yoyms Aq 
(061) ј95 јалопџеш 
әбе о әш L 


БЕНЕН ТЕС род —— Rl pus e a... i.s... sss =“ “ ==, — u ee — == 


| из | јовлејлороу JoBseL/LSs|soyuL 





Target turn generator. Switch threshold is set to 3 seconds. 


Figure A.12. 
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Figure A.13. Target and missile velocity computation. 
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Simplified 6DOF model with filter. 


Figure A.14. 
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Accelerometer. 


Figure A.15. 
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Figure A.16. IMU with additional outputs for tracking filter. 
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o—ÜB—Y tracking filter. Function block is ABGFILTER.M. 
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Full aerodynamic 6DOF model. 


Figure A.18. 
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Figure A.19. Aerodynamic moment and force models. Function blocks are 
ALPHABETA.M, AEROFORCES.M, and AEROMOMENTS.M. 
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Full aero model autopilot. 
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Figure A.20. 
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Full aero model 6DOF equations. Function blocks are EQNFORCE.M, 


Figure A.21. 


EQNMOMENT.M, EQNQUAT.M, and EQNPOS.M 
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Simplified 6DOF TBMD interceptor simulation. 
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Figure A.23. 
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Figure A.24. TBM thrust model. 
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TBM target model. Function blocks are BALLISTDYN.M and 


Figure A.25. 
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Figure A.26. TBM seeker model. 
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Figure A.27. TBM missile model switch. Function block is MODELSWITCH.M. 
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Figure A.28. TBM Missile Dynamics. Function block is GRAVITY2.M. 
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THIS PAGE INTENTIONALLY LEFT BLANK 


APPENDIX B. MATLAB@ CODE 


Filename Purpose 


|parasiticdragforce (|| 


draethesis parasitic drag force 


Table B.1. Matlab® Source Code Listing. 


Ji 







Filename 


w = a 


Table B.1. Matlab® Source Code Listing (continued) 
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function y=abgfilter(u) 
SABGFILTER Implements an alpha-beta-gamma filter as 


$ outlined in Bar-Shalom & Li "Estimation and 
% Tracking: 

% see also 

$ Copyright 1999-2000 by Triple B Enterprises 


кл» 
ххх хх 

%// File: abgfilter.m 

%// Мапе: LCDR Robert D. Broadston 

%$// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%/ / Compiler: MatLab v5. 3 

%// Date: 19 July 2000 


%// Description: Implements a 9-dimensSional state vector 

% / / alpha-beta-gamma tracking filter for use with 

% / / missile guidance laws requiring tracking filters 
% / / Note: uses global XLAST to preserve state 

% / / between iterations 

%// Imputs: measurements (los,los_dot,R,R_dot), 

%// missile pos (x,y,Z) 

%// Outputs: 9-dimensional estimate of target state 

%// теа v, wa z vz.az]"' 


%// Process:  alpha-beta-gamma filter outlined in Bar-Shalom & Li 


$// Assumptions: 
$// Warnings: may require up to 20 samples to stabilize from 


$// initialization 
БУ ӘӘК ЖСЖ ЖСЖ ЖЖ АСА ЖЭС ЖОТА ЖА Ж ЖТК ЖЖ АА КК Ж УЕ 


хххжхххххкхх 


%// Order of elements 


% / / -Define globals 

5// -Define constants 

$// -Define elements of input vector 
% / / -Functions 


БИЕ ЖАЖА ЖЖК ККЖ КАС ЖЖ ЖАС Жк КИК КА АРЫК ee 
ххххжххххх 


о Черт! >> 
global m d L XCG XCPN XCPW XCPB XHL 
global ST SW SPLAN SREF FILTSAMP XLAST 


$ ****** define constants ****** 


Gas ee define input- vector 555525 
JosSsdqot=1 (1) ; 

Pardoe=u (2) 

Jos h (53). 

phi=u(4); 

rdotzu(5); 

к= (Б): 

xm=u (7); 

ym=u (8) ; 

zmzu(9); 


25 


tee es Initialize variables са 
% compute target cartesian coordinates 
xt=R*cos (los) +xm 
yt=R*sin(los) +ym; 
zt=R*sin (phi) +zm; 


2 = Плус 2 21: 


% set noise parameters 
sigmav=1 ; 

Sigmaw=1; 
lamda-sigmav*FILTSAMP^2/sigmaw; 


$ set filter parameters from Bar-Shalom & Li 
Еајрћа=.9; 
fbeta-.9; 
Едапта=.9; 


% filter matrices 
-[1 FILTSAMP FILTSAMP^2/2 zeros(1,60); 


0 1 FILTSAMP zeros(1,6); 
0 0 1 zeros(1,6); 
2егоѕ (1,3) 1 FILTSAMP FILTSAMP^2/2 zeros(1,3); 
zeros(1,4) 1 FILTSAMP zeros(1,3); 
PN 5) 1 zeros (1, 3) 
zeros(1,6) 1 FILTSAMP FILTSAMP^2/2; 
Zeros (i) | FILTSAMP; 
zeros(1,8) axes 

= озо 0 о оо оно 
9:30:20 1520. 7050 DE 
Оооо 


$ compute steady state gains 
W=(falpha; fbeta/FILTSAMP; fgamma/ (2*FILTSAMP%2) ]; 


% build gain matrix 


P= [W Zeros (3,2) + 
Zeros (3,1) W zeros 1) 
zeros (3,2) W]; 


% * w * т ж Ж TF * Ж ж т Y+. Funct rons хх ххх 


новее 

xhet= bx Pues , 
хһасі-схһас-Р”(2-Нххһа); 
XLAST-xhatl; 

y-xhati; 


$//end o£ file abgfilter.m 
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function y-aeroforces(u) 
SAEROFORCES Computes aerodynamic forces on a missile. 


de de ge oe ge 


derived from Zarchan "Tactical and Strategic 
Missile Guidance" and Anderson "Fundamentals 
of Aerodynamics" 

see also AEROMOMENTS 

Copyright 1999-2000 by Triple B Enterprises 


K/L RK RR ЖКЖ ККК КККК КККК КККК ККЖ КЖК К ЖЖЖ АХ Ж 


* * * * * * * * * 


%// 
B 
%// 
n 
in 
su 
%// 
Е 
ЕЛ 
Su 
%// 
b Z 
Du 
%// 
%// 
E 


File:  aeroforces.m 
Name:  LCDR Robert D. Broadston 
MSEE/EE Thesis 
Operating Environment: Windows NT 4.0 Service Pack 5 
Compiler: MatLab v5.3 
Date: 7 Sept 2000 
Description: Computes aerodynamic forces for both subsonic 
and supersonic regimes on a symmetrical STT 
missile. 
Inputs: missile state, control deflections, angles of attack 
and rates 
Outputs: Body centered aerodynamic force components [Fx,Fy,Fz] ’ 
Process: Brute force computation of equations from Zarchan and 
Anderson 
Assumptions: 
Warnings: 


%//ХХХХХХХХЕХХХЖХХЖХХХХХХХЕХЖХХХКЖХХХХЕХХЕЖХХХХХЕЖХХЕЖХХЕХЖХХХЕЖХХЖжЖХХЖХАХЯ 


хжжхххжххх 


%// 
%// 
b 7 
Е 
%// 


Order of elements 
-Define globals 
-Define constants 
-Define elements of input vector 
-Functions 


Б/ / ЕТ КАКА А a a ова а и РОНА ГАР АЕ 


жхххххххх 


Би“ ~ + бартер ора е 5655s 
global m d L XCG XCPN XCPW XCBB XHL 
global ST SW SPLAN SREF 


$ ****** define constants ****** 


6 (oo Чар пе Input Е === 5 
States=u (L213) 

delta resu(l445 

delta e-zu(15); 

ЕЕЕ = Ш ds 

т атрпа=и (17); 

m beta-u(18); 

alphadotzu(19); 

betadotzu(20); 


о] 


% ****w** _yyyrtialize уаклташ ес pi ^ 


V om=sart (üu (2Wwo2+u(5)%“2+u (6) > $ missile velocity 
Machn=V m/machvalt (üu (3) ) ; % Mach number 
M_BETA=sqrt (Mach^2-1); % Beta factor 
(о ен о) ту m° 272; % dynamic pressure 
% хххххххххххх FUNCEI ONS ххх 
%----------- compute normal coefficients ------------------- 
%-------- these equations developed in Zarchan ------------- 
if (Mach>1.05) 
C Naz-2-3*SPLAN*m alpha/(2*SREF)... 
+8*SW/ (M_BETA*SREF)... 
+8*ST/ (M BETA#*SREF) ; 
(ние 2=8751/ (М ВЕТА"БОКЕЕ): 
C_Nz=C_Naz*m_alpha+C_Ndz*delta_e; 
C_Nby=2+3 *SPLAN*m_beta/ (2*SREF)... 
+8*SW/ (M_BETA*SREF)... 
+O*ST/ (M_BETA*SREF) ; 
C_Ndy=8*ST/ (M_BETA*SREF) ; 
C_Ny=C_Nby*m_beta+C_Ndy*delta_r; 
%------- these equations developed in Anderson ------------- 
else 
ов па |на: 
C_Ny=.5*m_beta; 
end 
%------------ compute drag -------------------- 
CD0scadO([states;thrust]); $ drag 
СВІ = Саз i NA Ny m alpha m беса, а Vm % coefficients 


drag= (CD1+CD0) *O* SREFE' ; 


%------------ compute forces ------------------ 
x=0; $thrust-drag; 

УС МО БЕ. 
ZC i NZ O~ SREE: 


М-ГЕ ЕУ ЕС |; 


%//епа ОБЕ file aeroforces.m 
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function y=aeromoments (u) 

SAEROMOMENTS Computes aerodynamic moments on a missile. 
derived from Zarchan "Tactical and Strategic 
Missile Guidance" and Anderson "Fundamentals 
of Aerodynamics" 

see also AEROFORCES 

Copyright 1999-2000 by Triple B Enterprises 


oe de de oe go 


"УДАЈЕ КЕ О ЕВО О О ВО АЈ ООСС КОЈЕ КО ЈОКО ЈОЈ 
* k k k kx x 


$// File: aeromoments.m 

%// Мапе: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler: MatLab v5.3 

%// Date: 7 Sept 2000 

%// Description: Computes aerodynamic moments for both subsonic 


$// and supersonic regimes on a symmetrical STT 

$ / / missile. Note: Moment about x-axis is 

% / / negative feedback of roll rate to stop missile 
% / / from Tolling. 


%// Inputs: missile state, control deflections, angles of attack 

% / / and rates 

$// Outputs: Body centered aerodynamic moments [Tx,Ty,Tz] ’ 

%// Process: Brute force computation of equations from Zarchan and 
DT Anderson 

%// Assumptions: 

%// Warnings: 


ЖКХ КАКА КК ЕЕК ООО ИЕС ИУ 
k kk kk xk ХХ 


%// Order of elements 


bu -Define globals 

5// -Define constants 

5// -Define elements of input vector 
%/ / -Functions 


Eff REXAKKKRARAAA AER RARE A ARERR ER RRR RAE BORN AAA RK ROR AK ARR A Ae аа ат 
cckock ck ko ko cko koX 

таста” че ше сора еи 

Global m d Lh XOG XCBNCXCPW ACPBOSUE 

global ST SW SEDAN SREN 


$ ****** define constants ****** 


о ее СПЕ vector r rinas 
states=u(1:13); 

delta г=ш( 4): 

delta_e=u(15); 

m_alpha=u(16) ; 

m_beta=u(17); 

alphadot=u (18); 

Depocot-u 0 9 
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БОРК irit t1alize varlas5l o sr 5 
V_m=sgrt(u(2)”23005)22+u (6) o) 

Mach-zV m/machvalt(u(3)); 

M BETA-sqrt (Mach^2-1); 

О Оха СПОТ ПС 272; 


missile velocity 
Mach number 
Beta factor 
dynamic pressure 


oP oP cO oh 


% хжхххххххххххх runctrions DX A хх ххх хл ххх 


%----------- compute moment coefficients -------------------- 
%------ these equations developed in Zarchan 
11 Масп>1.05 
C Myz2* (XCG-XCPN)/d... 
+3*SPLAN*m_alpha* (XCG-XCPB) / (2*SREF*d)... 
-9*5w* (XCG-XCPW)N (M _ВЕТЕЗЬВЕЕ С 
-8*ST* (XCG-XHL) / (M, BETA* SREF*d) ; 
C_Mdy=8*ST*(XCG-XHL,) / (M_ BETA*SREF*d) ; 
C_Mz=2* (XCG-XCPN) /d... 
+3*SPLAN*m_beta* (XCG-XCPB) /(2*SREF*d)... 
+O OW (XCG-XGPW) 7 (M_BETA*SREF ОСЕ 
+O *ST* (XCG=XHL) / (M BETA*SREF*dQ NX 
CUMGZ=5S ОТ). (т ЕРЕ СБВЕЕТО 


%------ these equations developed in Anderson 
else 

ОИЕ СЕЕ ОБ: 

ТЕ О 05: 
епа 


T x=-400*states (7); 


T y-(C. My*m alpha-«C Mdy*delta, e)*Q*SREF*d-800*alphadot; 
T z-(-C Mz*m beta-«C. Mdz*delta, r) *Q*SREF*d-«800*betadot; 


т 


%//end of file aeromoments.m 
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function y=alphabeta(u) 
SALPHABETA Computes angles of attack in both vertical 


% and horizontal planes 
$ see also 
% Copyright 1999-2000 by Triple B Enterprises 


Sf f k k k K k K X X X X k k X XK k % K £ K X Kk k X K KUK k UK X KUK X KO XK UK k K X X £ ххх хх Хх хх ххх ххх Ххх ххх жж 
KX KKK KK 

$// File: projxX.m 

$// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler: MatLab v5.3 

$// Date: 31 July 2000 

$// Description: Computes angles of attack using ATAN formulation 
$// in Bryson "Control of Spacecratt and Airerare: 
$// Inputs: missile state 

$// Outputs: angles of attack [alpha, beta] ’ 

$// Process: ATAN formulation of Bryson 

%// Assumptions: 

$// Warnings: 


qp f fco ke e ke ke ke e ke e e e e ce e eo e ОКО ОКК О ОВ КО А ЈЕ КОЈОЈ: 
*ockocko kck k kock k 


%// Order of elements 


$// -Define globals 

$// -Define constants 

$// -Define elements of input vector 
$// -Functions 


S f f kk kokck kk heo ok ke sk ke e ke ke ke ok ok oko ek oe kk ke ke ke koe ck X X X X < X X X X ke e e e kc kk ke kk kk Eck ck KC KK 
KKKKK KKK 


Е ше СО ае Е 

qoM exo define constants 595 5 

me tx xxx dese Еее "SuSE 
v= a ayuu e] 


$ ****** initialize variables ****** 


% < K Kk xk kk И Ae functions Kk xx xx xxx xx 

%---------- these equations developed in Bryson ------- 
$ using betal for sideslip angle to avoid problems with 
$ built-in matlab fxn beta 


alpramatctan2(v(3);sqx=ye (v (1) 21iv (2) ај 
betal=atan2(v(2),v(1)); 


y-[alpha;betal]; 


$//end of file alphabeta.m 
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BLL RRR KR KK KK ke he e e e KR e ehe che hee ke KK KK KK eR RK Ж Ж 


хххххх ko 


$// File:  auxplots.m 

%// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$// Date: 13 April 2000 

%// Description: Plots auxiliary variables from missile simulations 
%// Inputs: none 

%// Outputs: plots of auxiliary variables 

$// Process: попе 

%// Assumptions: попе 


$// Warnings: none 
Ф//ХХХХХЕХХХЕХЕХЕЖХХЕХХХХХЕТХХХХЕХХХЕХЕЕХХЕХХХХХХХХХХХХХХХХХХХХЕХХЕКХЕХХХХХҠЕЖЯЖ 


* * x x K xx 


%// Order of elements 


5/7 -Define globals 

É -Define constants 

БИМ -Define elements of input vector 
%// -Functions 


qo f f A kk ke hehe he ke hehe ke ehe See ke he hee ee he e ke e hehe e ek ke hehe ehe ehe echec hee he he hehe he heck ehe ce OX ORDRE 
*x * * x xxx * wx 


ЫС M deme globales tos 
$ ****** define constants  ****** 
psc el ne юре сог ~****= 


% ****** Initialize variables ****** 


% ** x xw wx ck ck ko ko functions ЖЖЖЖ Ж 


figure (4) 

subplot (4,2,1) 

plot (c Accelout/9:8045) 
ylabel('AccelOut') 
SsSubplot (4, 2, 2) 
plot(t,AlphaBeta) 
yYlabel(’AlphaBeta’ ) 
subpigt:Etd 23) 

plot (C, eúulers*57 3) 
ylabel(‘eulers’ ) 
sSsubplot(4,2;,4) 

plot (t,MissileV) 
ylabel('MissileV') 
ово О (4, 275) 
РТОЕ(Е, ССсетЕ-ЕОЕ) 
ylabel('AccelError’) 
subplot (4,2,6) 
plot(t,seeker) 
ylabel('seeker') 
Subplot(4, 2,7) 
plot(t,deltas) 
ylabel('deltas') 
$//end of file auxplots.m 


Lino tion = рапа (В) 
$B2QUAT Computes quaternions from a rotation matrix 


% 
% 
% 


B2QUAT (B) 
see also QUATERNION, BQUAT 
Copyright 1999-2000 by Triple B Enterprises 


qs f f Kok oko oe eoe e ee ee hehe he X X X ke e e e e ke hehe ee Se e e he e KX X XX XXX X X X X kX X X X X XX XX X e e e e e К 


ххххххххх 


BU 
Б 
%// 
%// 
%// 
%// 
БИ 
%// 
Би 
Е 
b 
ВИЙ 
%// 
%// 


File: b2quat.m 

Name: LCDR Robert D. Broadston 

MSEE/EE Thesis 

Operating Environment: Windows NT 4.0 Service Pack 5 

Compiler: MatLab v5.3 

Date: 12 Dec 1999 

Description: Computes quaternions from ABC rotation matrix 
using formulation of Kuiper "Quaternions and 
Rotation Sequences" 

Inputs: rotation matrix B 

Outputs: “quaternion аб ат а2 газ“ 

Process: Kuiper pp. 166-167 

Assumptions: 

Warnings: 


В Кю 


хххжххххххх 


%// 
%// 
%// 
БИ 
5772 


Order of elements 
-Define globals 
-Define constants 
-Define elements of input vector 
-Functions 


ок ны 


хххжжхххх 


О (Jerrie global s “ааа 


Фф ****** define constants ****** 


$ **^**^ define input Vector клк» 


$ ******  initialize variables  ****** 


% * * k * k * k *k xA X * * Funct onc хххххххххххх 


Сс^б=садагЕ( (1#В(1,1)+В(2,2у+В(3,3))/4); 


q lssqrt( 


1%В(1,1)- 2)=B(373)) 7/4); 


( B(2 
d 2=sgqgrt((1-B(1I, 1)+в(2, 2) BUS a ae, 
( ( 


q_3=saqrt ( 


(B 


a= 
b= 
c= 
d= ( 
e= 
f= 


( 
( 


J P (1 .1)=B(2 В 3)” O). 


в, Ы; 


ВВ, ВА: 
ОВ, 2) Виго шј) 4; 


ПНЕ А И НО, 


B( 
Biz, a) ЕЦ 2 d 
B ( 


Е 
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Dre acce De 


q 0--q 0; 
end 
if (a«0 & d«O 
qr e 1; 
епа 
if (b«O & d«O 
CC 
end 
1Е (с<0 6 е<0 
q 3=-c 3; 
end 


у= аа 


$//end of file b2quat.m 


& 


& 


с= 0) 


Е<0) 


е<0) 


Е<0) 
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Ешеісіе-лов у-ратсрі (14) 

%BANGPT Computes bang-bang control law for simplified 

% 6DOF model | 

% see also PROPNAVPT, VCPROPNAVPT, BRYSON, CHINGFANLIN 
% Copyright 1999-2000 by Triple B Enterprises 


ЛЕККЕ АСХАТ ККЕ КЖ X X X X K X X X X S Z kX X SASS ХК У ХИ 


* * *k k * ck ck ck ok 


%// Flle: BANGPT.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler: MatLab v5.3 

%// Date: 6 Aug 2000 

%// Description: Bang-bang control law for 6DOF flight model. 


% / / Uses two values of bang depending on range to 

E reduce problems with drag at start of engagement. 
5/7 -.005 rad/s dead band оп los rate 

%// Inputs: {seeker data,IMU data, timer] 

%// Outputs: {command accelerations,applied forces] 


%// Process: bang-bang control law 
%// Assumptions: попе 


$// Warnings: none 
Ак sos 


ххххххххх 


%// Order of elements 


$ / / -Define globals 

%// -Define constants 

у -Define elements of input vector 
Ф// -Functions 


Eft Rt RRR EKER AEA EAE EK REAR RR RR A KOR ARNE RR RRO Re Re ee 
ххххххххх 

и ее 9 ^^^ 

global m satflag 


ИК c derrne constants 757 
Nprimes5; 
Nprimez-5; 


ТЫ оне Тао ЕСО хе 
thetadot-zu(l); 
pusdot-ut20 
Tossu): 

ОАЕ а 
Е=и (6); 
Moms) 
heading=u (7); 
VIm= 0 (8) 
Wmdotzut9)- 
Dnhi=üu(10); 
Елпеса-ч ГІ); 
psi=u(12); 
сіше= и 1200 
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БИ А Ани" гла 1 те епа ес Ек» 
$ set max control force to 5 g (long range 
$ & 20g (end game) 


if (R»5000) 
Nbangz5*9.8045; 
else 
Nbangs20*59.8045- 
end 


9; Kk oko ck x x ko ko ko tunctrone хххххххххххх 
$ establish a dead band on theta dot 
if (abs(thetadot)».01*pi/180) 
ny-Nbang*sign(Vc*thetadot); 
else 
пу=0; 
епа 


nz=Nprimez*Vc* (phidot) -9.8045; 


% control force limiter 
ЈЕ satflag 
а ере (за == 0“9' 8045) 
ny=signi (ny) "3079.8045; 
end 
if (abDs(nz)=309.8045) 
па = ено та зао о а Бе 
епа 
епа 


% compute ABC forces applied 
Jape 658 

Fy=ny*m; 

Fz-nz*m; 


% output vector 
y= lay: nz Bx Баве |; 


$//end of file BANGPT.m 


106 


function y=bryson(u) 
S$BRYSON Computes optimal guidance law derived by Bryson & Ho 


% with dragforce inputs for point mass simulation 
% see also PROPNAVPT BANGPT CHINGFANLIN 
% Copyright 1999-2000 by Triple B Enterprises 


% / / * * * * k k k k * k k * k k k * k * * + k Kk kk kk kk 'k'k'k kk kk kk kk kk k kk k k k * k kk kk kk oe 


ххххххххх 


$%// File: »vson.m 

$// Name:  LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windews Ni 4,0 Service Раск = 
%// Compiler: MatLab v5.3 

%// Date: 18 Sep 00 

%// Description: Modified PN differential games guidance 


%// law from Bryson 6 Но 
%// Inputs: Seeker outputs, filter outputs, missile timer 
%// accelerometer output 


$// Outputs: command accelerations, y and z forces for drag model 
%// Process: 

$// Assumptions: 

$// Warnings: 


% / / * ** k 'kk'Kk'k'k'k k kk kk k kk k kk kk kk kk kk kk k * * * * e ek hee e e kc se ce e e ke o e ke ke e c ke kc ke ec ke 
ххххххххх 


2// Order of elements 


%// -Define globals 

$// -Define constants 

i -Define elements of input vector 
%// -Functions 


$ / f koe koe e e ck ehe e e e Se ke ke ehe ee Kk e ke ck cce e e cce e e e ke ce ke e e Se he ce ce e ck ck ck che E ke koe kk e e kk ke ke kk 
* * k * K * * ЭЖ 

Б бепе дора кке, 

global m satflag 


Б "*ocdertine constaHts 15 E 
Nprime=3; 
Nprimezz3; 


И - <= <=“ define input vector === 
thetadot=u(1); 
phidot=u(2); 
los=u(3); 
рое 
R-u(6); 
ен) - 
heading=u(7); 
Vis). 
М\тсос=п (О); 
phi 10). 
thetazu(11); 
рее): 


п =Ссасе= 15:21); 
timezu(22); 


асолати 207725); 
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ФО Рики initialize variables ****x* 


% x * xx xx x x Хх ~ x une So ms ck ck ck ook oo k ñ kk 
1f time<2.0 
ny=Nprime*Vc* (thetadot) /cos(psi-los) ; | 
nz=Nprimez*Vc* (phidot)/cos(theta-philos)-9.8045; 
else 
eGp=50*9 78045; 
ce_lat=sqrt(m_state(3)%*2+m_state(6)%2); 


if (ce lat==cp) 
се 1ас=29 ^9 8045: 
епа 


ce vert-abs(m state(9)); 


if (ce vert=-Cp) 
ce _ vert=29*9 28045; 
end 


ny3 (l ce lat/cp)*Ve*thetadot: 
nz=3/(l-ce vert en) Ve'phidot-9.8045; 
end 


if satflag 
if (abs (ny)>3079. 8045) 
Е Оо. 
епа 
if (aDs(nz)>30*9.8045) 
nzssign (nz) о а: 
end 
end 


Fx=0; 
Fy-zny*m; 
Tuc cm mc 


ра 


senc oS 1те ркугоп m 
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Luncetuonsvedoiu 
$ CDO Computes induced drag coefficient 


% 
% see also CDI 
% Copyright 1999-2000 by Triple B Enterprises 


gs f f| Noc ok ck oe e oe eoe oe oe e ehe e ce e e e e e e che e ke ce e ke e ce ce e e cce ee kk e Sec kk kk kk k e koe 


*oc*ockockock k k k Ж 


лете: сабт 

%// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

$// Date: 9 May 00 

%// Description: computes induced drag coefficient for full 
$// aero model 

$// Inputs: state, boost status 

%// Outputs: drag coeffient 

$// Process: polynomial fit to data from Hutchins EC4330 notes 
%// Assumptions: 


%// Warnings: 
Q; f f koe ke ehe ke oe e e e e ke e e e ce e e e e e he he he ke e e e ck e e e See e e e oe he he ke ee e e ke ke e Se e ke e khe ee e ke e ce ke ke e ko 


ck ck ko ck k kk 


%// Order of elements 


% / / -Define globals 

%// -Define constants 

%// -Define elements of input vector 
5// -Functions 


Г ROCK CK CKOK X KX X X KX K X X X X X X X K X X XX KX X Ж ж Ж X X X X X X X X X X ES аа аа атала а а 
ххххххххх 


ЫК к=к» ре о mran 


ф ****** define constants ****** 
NoBoost=[-0.0014 0.0299 202900 0.26256]; 
Boost=[-0.0012 0.0243 -0.1521 0 4044]; 


~ “~ -~> сеттпе snput хебког |“ | = 
ЕСЕ (4) 2+0 (5) 2+006) 2); 
alt=u(3); 


boostsutld)- 


$ ******  Jinjtialize variables ****** 
mach-v/machvalt(alt); 


% хххххххххххх РИПСЕТОПЕ * *& k k k Kk x= ko ko ko ko 


ІТ (шасп>100) 
mace. 3 
end 


109 


% these curves approximated from "typical" data presented 
% in Stevens & Lewis and Hutchins 


% compute Cd0 

if (boost & (mach<1) ) 
y=.15; 

end 


if (~boost & (mach<1) ) 
у=.25; 
епа 


lif ((mach>=]) & (boost--0)) 
y-polyval(Boost,mach); 
end 


if ((mach»-1) & (boost--0)) 
y-polyval(NoBoost,mach); 
end 


if ((mach»5) & boost) 
y=. 10; 

end 

if ((mach>6.4) & ~boost) 
УЕ Во. 

епа 


$//end of file cd0.m 
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Пао у= са (1) 


$CDI Computes induced drag coefficient 
% see also CDO 
% Copyright 1999-2000 by Triple B Enterprises 


2 / / X XK X X X X X K K k K X X K X RK ххх e ke e e ke ke ke ke I I Хх хх хх хх ж хк 
жжхжххххххх 

Е сам 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 18 Арк 00 

$// Description: computed induced drag coefficient for full 
ъ// aero model 

%// Inputs: see below 

р ОЕ" cdi 

%// Process: Anderson 

%// Assumptions: 

%// Warnings: 


$/ [ AA AAEE EKK ESE K A e e e e e X X X X X e e e e X X XX X X ke e he e e he e ke e sk e ke КЕКЕК ККК ККУ KES 


ххжхххххх 


%// Order of elements 


%// -Define globals 

6/ / -Define constants 

%// -Define elements of input vector 
%// -Functions 


Qo f ooo e e ehe e X X X X X % he ke e e ke e X XIX X X X hehe X e he ke coke e e e e 2 X ke e e ee ke koe e ck ke he e ke cec e e e e e e e 
ook oko xk xx x x 


Бе Е: пе орайы жя 
$ ****** define constants ****** 


и efine Input vector 57552 
ч л=п (у; 

C_Ny=u(2) ; 

m_alpha=u (3); 

m_beta=u (4); 

elles (sie 

v=u (6); 

Ser initialize variables r ri nn 
Mach=v/machvalt (alt); 


$ *ockock ck coc X x A X kx xk ЕпревјОЛе жххххххххххх 
%-------- these equations developed from Anderson ------- 
if (Mach>1.0) 
M BETA-sqrt (Mach^2-1); 
Cdi-s(4*m alpha^2/M, BETA«4*m beta^2/M BETA); 
else 
Cda= (Cony? 246 NZ eer 
end 


Ба 


%//end of file cdi.m 
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function y=cdvmach (mach, boost) 
$ CDVMACH Computes approximation of zero lift drag 


% coefficient vs. mach number 

% CDVMACH (MACH, BOOST) 

% see also MACHVALT 

% Copyright 1999-2000 by Triple B Enterprises 


Bf [ERE RRRRRAEKRKREKREREKRK EKER RK KE KE EA KK KK eee eee 
KKKKK KK KK 


$// File: cdvmach.m 

%// Name:  LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$?// Date: 17 Арг 00 

$// Description: computes polynomial fit for cd0 vs Mach number 
%// Imputs: mach # and boost status 

ФІ Outputs: сао 

$// Process: Fit on data from Hutchins EC4330 notes 

%// Assumptions: 


%$// Warnings: 
о о оно 


ЖЕКЕ ЕККУ 


%// Order of elements 


%// -Define globals 

% / / -Define constants 

5// -Define elements of input vector 
%/ / -Functions 


KL fRRERKRKKRKRERKKEEKEKEREKRKEREKKKEKKEKKEKKKKERERREKKEKKKKKEKKKK KKK KKK XX ы 
* * k K k * KKK 


о р о tM 


Еее ЕЕ 5 
NoBoost=[-0.0014 0.0299 -0.2110 0-6256]: 
Boost=[-0.0012 0.0243 -0.1521 0 4044]. 


ыыы ыды JO Fir O Uy ооо C Or 2722 


ф ******  initialize variables  ****** 


% * * k £ k k k k & k * Ж fumetronpDs хххххххххххх 
її (роозе & (mach<I)) 

= >; 
епа 


1£ (~boost & (mach<1) ) 
i AE 
end 


1f ((mach>=1) & (boost~=0) ) 
y=polyval (Boost,mach) ; 
end 


if ((mach»-1) & (boostzz0)) 
y=polyval (NoBoost, mach) ; 
end 


Шы ((Пасп>о) 5 роовб) 
у=. 10; 
епа 


if ((mach»6.4) & ~boost) 
у=- 132; 


епа 


%//end of file cdvmach.m 


funGcEren y- Chingfanl iny 
$CHINGFANLIN Computes optimal guidance law derived by Ching Fan Lin pg. 
ЕШТЕ 


% with dragforce inputs for point mass simulation 
% see also EXACTPROPNAV2 
% Copyright 1999-2000 by Triple B Enterprises 


qs f f Kok koe e e e he e eee e e ee he e e e e he ee e e e ce he e e e e e e Se ehe e Se e e ke eee he ee ee e e Se ce ccc ck e e heh kk 
ххжххххххх 


ое. сїїпдаїап11п.ш 

$// Маше: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler:  MatLab v5.3 

%// Date: 15 Sep 00 

$// Description: computes APN guidance law from Ching Fan Lin 
%// Inputs: seeker output, filter output, accelerometer, 

% / / missile timer 

$// Outputs: command accelerations, y and z forces for drag 
%// Process: 

$// Assumptions: 

$// Warnings: 


$ / f[ Kok e e eoe e e e ke ee x he he k K K e e e ee ee e ee e ehe e e ee СКС ССК ke cce ke ck c e e ce ce e e kc koe 


жхххххххх 


%// Order of elements 


$// -Define globals 

% / / -Define constants 

$// -Define elements of input vector 
%// -Functions 


% / / * * * * * * * k k k k k e K k k R KOK e kO k ko ko k ko kO k K k ko K e e he sek KO ko kok R ko k Sk Ko K KO hee hh Se Se Se sees RR RO SO OO 
ххххххххх 


$ ее е ях 
global m satflag 

S **5*** ЕТ СЕНЕ "жс 
Nprime=3; 

Nprimez=3; 


$ тихи Чер тешр әр шү EX XX 
певао в= а s 
phidor=- u car. 
lossut3)e 
РО 
R=u (6) ; 

Ve=-u (5); 
heading=u(7); 
Vm=u (8) ; 
Vmeot=u 9): 

Ко = АЕ: 
theta=u(11); 
Р-НЕ), 


tgt_state=u(13:21); 
time=u (22) ; 


aceel тп=п(23:25):; 
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% ****** initialize variables ****** 


% Kak kkk Kk e e xx functions ххх Хаж 
ШТІ (Ус--0) 
tgoz1e6; 
else 
EGO=n7/ Vc; 
end 
$ compute relative state estimate 
хпас=[Е*со5 (16656); 
F onn (Tos) 
БУБА рО 06) 
tgt state(2)-Vm*cos(psi); 
tgt. state(5)-Vm*sin(psi); 
tgt state(8)-Vm*sin(theta); 
tgt state(3); 
tgt state(6); 
tgt_state(9); 
accel inii); 
ассе1_1п(2); 
accel_in(3)]; 


if time<2.0 
ny=Nprime*Vc* (thetadot)/cos(heading-los); 
nz=Nprimez*Vc* (phidot) -9.8045; 


else 
ucz(5/tgo^2)*[eye(3),tgo*eye(3),tgo^2/2*eye(3),zeros(3)] *xhat; 
nysuci2); 


mz=ue(3)-9-6045- 
end 


if satflag 
if (abs (ny) >30*9.8045) 
ny=sign (ny) *30*9.8045; 
end 
lf (арабп2)>30%9 5045) 
па=елаи па) 3075. 5045: 
епа 
епа 


Ех=0; 
Fy=ny*m; 
EZ=AZ mM: 


ү- (пу;пт””ЕХЕМ;Е2І; 


$//end of file chingfanlin.m 


ІШЕ 


function y=draginduced(u) 
SDRAGINDUCED Computes induced aerodynamic drag force 


% 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


Bl [RRR RK RK KR KR RR RR KR KK KR KK KK RR KK KK KR RR KR KR KK RK RK KK KK KK RK KR KK KK KK KK KKK KK KK 
ххххххххх 

$// File:  draginduced.m 

%// Мапе:  LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler: MatLab v5.3 


$// Date: 19 Sep 00 

%// Description: computes induced drag for simplified 6DOF 
%// Inputs: force output of guidance law, state 

$// Outputs: drag force 

%// Process: work backwards to CN from forces 

%// Assumptions: 

%// Warnings: 


BLL RRR RE RERKREEKEKRKEKREKRKKEKREEKRKEEKKREEEKEKERKKRKEKKKEREKKEREEKKEKKEEREK KK KAR 
ххххххххх 


2// Order of elements 


$// -Define globals 

%// -Define constants 

577 -Define elements of input vector 
$// -Functions 


S f Kok kk A kk ok X X K X K K. K K UX KOK X X KUK X X X X KS E S X X K X k AX X K k X X X X k K k X X X X K Z 2 2 XA & Kee 
ххх 

5 кек -детпе" дадшорај м 

global SREF m 


$ ****** define constants ****** 
еАК=1.5; 


oe 


elliptical eff & AR 


o xa ыар т түбө тр лесе 2222: 


Fy=u (2); % у force 

Fzzu(3); $ z force 
v2suo2)92v*0(5)^2:0(9)^2-5 $ missile velocity 
alt=u(6); $ missile alt 


ЕЗ ә "паса се сакавте e Et 


rho-zrhovalt(abs(alt)); $ atmospheric density 

mach-sqrt(v2)/machvalt (alt); 

(= потекла 2; % dynamic pressure 
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% жххххххххххх fonctions хххжхххххххххх 


if (Q==0) 
Cny=0; 
а= 07: 


else 
CnyzFy/(Q*SREF); $ y normal coefficient 


Фоа Јел А СО а % z normal coefficient 
end 


co 


Cdi= (Cny22+Cnz22) 7 (Di *expy induced drag coefficient 


if (mach<1) % subsonic drag ecual to 
Cdi=.25*sqrt (Fy*2+F2%2)/(m*9.8045); % max CdO*applied G force 

end 

pu cS ESTIS $ drag force 


$//end of file draginduced.m 


ТІ? 


function y=draginducedtbm(u) 

SDRAGINDUCEDTBM Computes induced aerodynamic drag force 
% 

% see also 

% Copyright 1999-2000 by Triple B Enterprises 


Q5 f [KK ke ke ke ke ke ke he ee ke ke КЕ КЕ КЕ 
* x x s k Xe 

$// File: draginducedtbm.m 

%// Мапе: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 


$// Date: 19 Sep 00 
$// Description: computes induced drag for TBMD simulation 


$// Inputs: force output of guidance law, state 

$// Outputs: drag force 

$// Process: work backwards to CN from forces, corrects for ECI 
%// Assumptions: 

$// Warnings: 


% / / * ** oe e ke ke ke ke eoe e e Se ee ke ke ke ke ecce ke ck ke e ecce ke e ke e ce ke e ce ck e ck e e e ke ke ke e ck ke e ke ke ke e ke e ke ke ke e kk kk 
oc * < * AX Ax * x 


%// Order of elements 


% / / -Define globals 

$ / / -Define constants 

% / / -Define elements of input vector 
% / / -Functions 


% / / eoe e e e ke e ke ke ke ce e ke e ke eee ce ee e ke ee ke ce ck ke ke ke e ke ecce ke ke e ck e ce ck ke e ck kk ke ke ke e ke ke ke ke ke e o ck ck ke kk ХХ Ж 
* * *K X ck ck k kc 

Тат сүс шешенке ыс o “еқ 

global SREF m 


$ ****** define constants ****** 
eARz1.5; 


со 


elliptical eff & AR 


АНА define Ine. vector "яи 


Fyzu(2); $ у force 

Ez-u 30 % z force 

МАЕ Тан о). $ missile velocity 
alt=sqrt(u(1)%2+u(2)%*2+u(3)%2)-6371e3; % missile alt 


$ ****** Initialize variables ****** 


rho=rhovalt(abs(alt)),; $ atmospheric density 
mach=sqrt(v2)/machvalt (alt) ; 
О= патот $ dynamic pressure 
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9; oko ck xe Ax x x kx x x X functions ockock ck ck ck k k *< K < *< 


1Е (О--0) 
Спу-0; 
Cnz=0; 
else 
Спу= Рума * 5КЕЕ):; ишао сте Е ое ВЕЕ ТЕЛЕ 
(прва (Or SREF) % z normal coefficient 
end 
Сат=(Спу^2+#Сп2^2) / (pi 'eARBR)5 $ induced drag coefficient 
if (mach«1) $ subsonic drag equal to 


Сал= 25*&sgrt(Fy^2t-*Bz^2)/ (m9 528045); + па“ саб“арр лева С ЕОс 
end 


y-Cdi*Q*SREF; $ drag force 


$//end of file draginducedtbm.m 


ll. 


function y=dragthesis(u) 
S$DRAGTHESIS Computes aerodynamic drag force 


% 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


Qf [Ree ke eoe eee ce ke ke ke he e cce ccce ke ХХ» ххх» ХХ ххх хх» Ж ххх» ХКК ЖХХ Ж» ЖЖ ххх Ж ЖЖ К К У 
ххххххххх 


$// File: dragthesis.m 

%// Мапе: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$// Date: 19 Sep 00 

$// Description:  computes parasitic drag after breaking apart 
$// state vector 

$// Inputs: state vector, boost status 

$// Outputs: parasitic drag force 

$// Process: 

$// Assumptions: 


$// Warnings: 
$/ / Ck oe ce ee e ke Sk ke ke he he he ke e ce ck ke ke e cec ke ke ke e e ke hee e ce e cce ck ck ck ce ce ce e ce ce ck ck ck ck cock ck ck kk ck kc ck kk ck ko kk ko 


яхххххххх 


%// Order of elements 


%// -Define globals 

% / / -Define constants 

% / / -Define elements of input vector 
% / / -Functions 


85 f [Ck ke ke ce e ke he he e he hehe ke e he ke ke ke e e ke ke ek ЖАКА АЖАН А УУ ЖАКА ЖУАН А АЖАН ХК 


хас хх ж 
p rre Ее бора ===> 
global SREF 


БРАШНА и ое 75-5 
domes deblsucelupub vector “о 
уе12=0 (4) ^2+0 (5) ^2+0(6) ^2; 

al uo): 

boost=u (14); 


$ ****** Initialize variables  ****** 


% хххххххххххх те топе ЖАС Ж ЖТА ЖС 


y=formdrag(SREF,alt,vel2,boost); 


$//end of file dragthesis.m 


gb f kk хх кокк кк хх к хх хх кх ххх хх клх КК к КЕ Хх e ese cce e eee hee cce Ж 
ххххххххх 


$// File:  drawmissile.m 

$// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

$// Date: 31 May 00 

%// Description: draws picture of current missile defined 
$ / / by missiledata#.dat 

$// Inputs: 

$// Outputs: 

$// Process: 

%// Assumptions: 


%// Warnings: 
EE A EAA AA E A E A E A К 


Ck ck * ko ck ck ck o 


%// Order of elements 


%// ~Define globals 

%// -Define constants 

$ / / -Define elements of input vector 
$ / / -Functions 


Sif м ххх КО АЕК ххх ххх кк кх ЕКА Ххх хк КККК ЫККА КА КЫК Аль 
Ok oko ck ck koX x 


зея дељапе сјорваје “ык 


$ ****** define constants  ****** 
ire Мола 25-013 
cepe us. 75 5. 7512 


ле ше ое ~ es 


5 “ХЕ лл латте variables ass s. 
позех=[0 LN LN 0]; 
повеу-(0 d/2 -а/2 01; 


Poco [PN L L EN LN]; 
ову у- [472 а/2 -а/2 -а/2 9/21, 


wingx-[LN-XW LN+XW+CRW LN+XW+CRW LN+XW+CRW-CTW LN+XW LN+XW] ; 
wingy-[d/2 d/2 d/24WXT«HW d/24WXT-HW d/24WXT d/2]; 


tailx-[L-CRT L L L-CTT L-CRT L-CRTI- 
taily=[d/2 d/2 d/2+TXT+HT d/2«TXT«HT d/2+TXT d/2]; 


CPEFF=(XCPN*AN+XCPB*AB+XC PW* SWtXHL*ST) / (AN+AB+SWtST) ; 


$ compute time 

time=rem(now,1); 

пр вроре Бете“ 24). 
Mins=floor(rem(time*24, 1) *60) ; 

Gimestr={=(° ~,mumZstri ht)” > enum БЕ Е ое 


% хжхххххххххх Fun E teo ms хххххххххххх 


figure(10) 

eT 

пота on 

axis equal 

fill (nosex,nosey, ‘w’) 

Емба bodyy, magne.) 

fill (wingx,wingy, dgrey) 

fill (wingx, -wingy,dgrey) 
fill(tailx,taily,dgrey) 
fill(tailx,-taily,dgrey) 

РОВ 0, КО’) 

plotoshm 9 r*') 

DIOC ОРЕКЕ О, а“) 

legend(‘Center of Gravity’, ’Hinge Line’,’Effective Center of Pressure’) 
title(['Missile Plan View ',date,timestr]) 
xlabel('meters') 

ylabel(’meters’ ) 

hold off 


%//end of file drawmissile.m 


function y=dynamic3D(u) 
SDYNAMIC3D Computes motion dymanics for a 


% body in three dimensions 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


Bl [RRR KRKKKKKKKEKREKKEKKEKEKKKEKEKRKKRKKKKKKEKEKKKKKKKEKKKKKKKKKKKKKKK KKK KKK KKK KKK 


KaKKKK KK KK 

$// File: dynamic3d.m 

2// Name: LCDR Robert р. Broadston 
%// MSEE/EE Thesis 


%// Operating Environment: Windows NT 4.0 Service Pack 5 


%// Compiler: MatLab v5.3 
%// Date: 16 Feb 00 
$// Description: target motion dynamics 


$// Inputs: target state, turn rate input 


%// Outputs: derivative of target state 
%// Process: 

$// Assumptions: 

%// Warnings: 


Bf BAAR RAR RK RRR KN ERR RK RRR RRR KK ERR KAA BORK AK A AK NR RRR ARR RR ла 


* * x w w KKK 


%// Order of elements 


%// -Define globals 

$// -Define constants 

5/7 -Define elements of input vector 
$// -Functions 


а о а о о о о ооо ое ое 5 


ххххххххх 


Би“ === ~ „де тпе. доба о кон ` = 
$ ****** define constants ****** 


osek define input vectores osse 
omega=u(1); 

о»: 

и): 

y=u (4); 

ео. 

zm 5 

Ер; 


$ ****** initialize variables ****** 


% хххххххххххх functions хы k * K K < x x x 
peo: 

-omega*ydot; 

VIOE; 

omega*xdot 

ОЕ 

О 


%//end of file dynamic3d.m 


function y=eqnforce(u) 
SEQNFORCE Computes force dymanics for six degrees 


% of freedom flat earth model 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


qb f [eee ck ck Sk e e ke e ke ke ke ke e ce e ee e e ce ke KR KKK KKK KK RK KKK KK KR RK KK KKK KK KK KK 
oc ck ck ok ko ok ok ox 


$// File: eqnforce.m 

$// Маше: LCDR Robert р. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

%// Date: 10 May 00 

$// Description: force dynamics for full aero model 
%// Inputs: see below 

$// Outputs: solution to force equation 

%// Process: 

%// Assumptions: 

$// Warnings: 


%/ f e ke e ke e ke ke e ke ke e ke ke ke ke ke ke ke ke e ke k k k K k K kk kk k kk kk k kk kk kk kk kk kk kk kk kk kx 
kak kkk KKK 


%// Order of elements 


% / / -Define globals 

$ // -Define constants 

$// -Define elements of input vector 
Ф// -Functions 


Sb Jf ok Doo ke ke ke ee e ee Kk * ke ke ke ke KR ke ke e e ke ke ke ke ke ce ke KR KK KK RK KK RK KK KR KK KEK KKK KKK KKK KR KK KKK 
ххххххххх 


В ТТ 5 
$ хо пене сопевелп еы ^к» 


те vector = 
МЕ ООО ОЗИ 

ЕЕЕ ООо 

omega b= un s o n o); 

P=u(7); Q=u(8); R=u(9); 


а= ОСО С u looo 


пада=ѕагі (а (1) ^2+9(2) ^2+4(3) ^2+9(4) ^2); 
а=а/ таса; 


% the ever lovin’ force of gravity 


а лы КЕ 
% note we are not using an external gravity model here 


= РО РО ОБЕ Е 


% and lest we forget, mass 
m=u(17); 


ок initialize variables ****** 


* * k kX kk kk k xk k 


$ хжхххххххххх functions 
$ some heavy duty number crunching 
$ compute rotation matrices 
Bequat 2 od): 
OMEGA_B=[0 -R Q; 

R O -P; 

SON PIE: 


y=-1*OMEGA_B*v_b+B*g+(1/m) *F_B; 


$//end of file eqnforce.m 


function y=eqnmoment (u) 
SEQNMOMENT Computes moment dymanics for six degrees 


% of freedom flat earth model 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


Kf [RRKRKRKKKEKEKRKEKKKEEKKEEKEKEKKKEREKKEKEKEKEKRKEKKEKEKEKEKEKKKKKKKEKEKE KK KERR KK EX 


хжхжхххжххххх 

$// File: eqnmoment.m 

%// Мапе:  LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$// Date: 11 Sep 00 

%// Description: computes moment dynamics for full aero model 
%// Inputs: see below 

%// Outputs: see below 

%// Process: 

%// Assumptions: 

%// Warnings: 


ААА Ак ЛК К eee 


ххххххххх 


%// Order of elements 


% / / -Define globals 

%// -Define constants 

% / / -Define elements of input vector 
%// -Functions 


MX XX X XK ао AK ok а аа тата ты 


ххххххххх 

ккк» чар пе glopals SES 

qox Ex define constants ***xxx 

$ wx xx define іпрНнс voc Or "= +.“ 
omega_b=([u(1);u(2);u(3)]; 

P=u(1); Q=u(2); R=u(3); 


% torques 
ТАВ- а Е Е: 


$ inertial matrix 

J=[ü (7), 050; 
(лане 
Ung CS - 


% апа lest we forget, mass 
mzu(10); 


6 КАНАЛА JDT а Га гема ево ПЕРА АН 
$ ххх ххх xXx k functions хххххххххххх 
% some heavy duty number crunching 
OMEGA_B=[0 -R Q; 

R 0 -Р; 

ТІ ШЕРДІ 


y=-1*inv(J) *OMEGA_B*J*omega_b+inv(J) *T_B; 


$5//end of file eqnmoment.m 


function yzeqmposit (u) 
%EQNPOSIT Computes NED dymanics for six degrees 


% of freedom flat earth model 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


Bf [RRR KRKEKKKKKKEKEKKKKKKKKKKKEARRAKKKEKAKKKKEKKKKKKEK KA KKKKEKKKKKKKKK KKK KKK * 
kek kkk kkk 

$// File: eqnposit.m 

%// Мапе:  LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

$// Date: 10 Мау 00 

%// Description: navigation equation for full aero model 
%// Inputs: see below 

$// Outputs: inertial velocities 

$// Process: 

%// Assumptions: 

$// Warnings: 


% / / > k kk Kk kk kk kk K kk kk k Kk kk kk ok e ee ke sk e e sese ee se oec oko ek kk Gk ko Rock Ж ЖЭС Э ЭЖ 
хжххххххх 


$// Order of elements 


%// -Define globals 

%// -Define constants 

$/ / -Define elements of input vector 
$ / / -Functions 


Q; f f[ k * * * * * * k k W k k k k k k k kO kO SK e esee eee ke e kO Sk kO kO ko ke ok ee ee e e k KO ko Kk kO e Sk e e e Ko Ko Ko koko kO kO Ж KO KO KO O O K 
* * * + ХХХ 


ВА ен derine glopals gme А а 
$ ****** define constants  ****** 


р derine пр ес ог ш кук 
В: 


ае а Сау Ој па по јела (0 


% ****w** Initialize variables ****** 
Magd=sqre (ail) “240(2) “24a: 259040) 2 
q-q/magq; 


% хжххххххххххх functions хжжххххххххх 


% compute rotation matrices 
B=quat2b(q); 


у=Вв “у р; 
$//end of file eqnposit.m 


function y=eqnquat (u) 
$EQNQUAT Computes quaternion dymanics for six degrees 


% of freedom flat earth model 
% see also 
$ Copyright 1999-2000 by Triple B Enterprises 


ккал ee 
KK KKK Kk KK 

%// File: eqnquat.m 

%// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

$// Date: 13 Sep 00 

$// Description: computes quaternion dynamics for full aero model 
%// Inputs: see below 

$. Qülputs: “wq dot 

%// Process: 

%// Assumptions: 

%// Warnings: 


Е Е К Аа гаса иа 
хжхххххххх 


%// Order of elements 


% / / -Define globals 

%// -Define constants 

$ / / -Define elements of input vector 
$// —Functions 


QJ Kk X K k Ok K k k k k UK K X K K K K K K k k k kk K k k KOK OK K KOR KOK K KOK SSK k КОКК К КОК ААА КОК КОККЕ ОКЕ САГА ААДА 
o oc A x kx kx x 


5 а Е. НЕ ОЕ 55205 
фени Cometanta 4 ***xx 

Go кеке ОЧСО Б INPC Vector "rinra 
с-ш “а атак Ра АТ; 


omega_b=[u(5);u(6) 


am zu ie 
В ОЕ -uuE 


(57) 
aD 
9 ххх „ара -2е variables  ****** 
magq=sqrt(q(1)*2+q(2)*2+q(3) *2+q(4) %2); 


q-q/magq; 
% ххххххххххкхкх functions kek KK kkk Kk KK 
OMEGA q-[0 P О в; 

-P 0 -R Q; 

-Q R 0 -Р; 

ӘН Oo CPI 


qz-(1/2) *OMEGA q*q; 
maggesemctgui)s2-qu2)99 ceo) 
if (magq--0) 

y-q/magq; 
else 

УЕ, ОО: 
епа 


$//end of file eqnquat.m 
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function y-flatearthdyn(u) 
$SFLATEARTHDYN Computes motion dymanics for six degrees 


% of freedom for a flat earth model 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


Жк кк кАк ххх ec CK koe ck e ok ek ok ke А хх Жк ke хк 
ххххххххх 


%// File: flatearthdyn.m 

%// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

$// Date: 1 Aug 00 

$// Description: computes 6DOF dynamics for flat earth using 
$// quaternion formulation 

%// Inputs: see below 

$// Outputs: derivative of state vector 

%// Process: Stevens & Lewis 

%// Assumptions: 

%// Warnings: 


а ee 
kkkkkkk kk 


%// Order of elements 


$// -Define globals 

$// -Define constants 

%// -Define elements of input vector 
$// -Functions 


eh fae RRR ARK RRR RN A eK Xe KKK RR ER RAR RRR ARK KR RR RK KN i 
* * * x Kk kK KK 


au ce" define слорашенн “=== 
Z ****** define constants ****** 


X ***x*- define Inpa Vector <***** 
p [u (1); u 2) ale: 

ЕЕ Оа) Оооо ИИ 

omega, bz[u(7);u(8);u(9)]; 

P-u(7); Q=u(8); R=u(9); 


а= о ооа зт: 


пааа= Зате) 2а 2 е(з) 2-а (4) 2). 
q-q/magq; 


>= p: о оте са nd]. 


% inertial matrix 

= паша 05 о 
БАП 
ОО 6 


$ forces 
БЕТІ (ЛӨ) utto] 


$ torques 
Pee (20) eu 2) (220 


% the ever lovin’ force of gravity 
$ note we are not using an external gravity model here 


а=[0;0;9.8045]; 


% and lest we forget, mass 
m=u (26); 


+ “к=к a bere emote Ores. Хх 
$ compute rotation matrices 


В=еџа 20 (а); 
OMEGA B-[0 -R О; 
R 0 -Р; 
-O P 0]; 
OMEGA q-[0 P Q R; 
-Р O0 -R 0; 
-0 R 0 -Р; 
-R -0 P 0]; 
% хххххххххххх functions хххххххххххх 
у= [ zeros(3), В”, zeros (3), тегов (ЗЕ 
zeros(3), -OMEGA B, тегов 5), zeros (3748 
zeros(3), zeros (3), -1*inv(J)*OMEGA_B*J, zeros(3,4); 
zeros (4,3), zeros(4,3), zeros (4,3), (-1/2) *OMEGA q]; 


у=у*х; 

y=y+[zeros(3,1); 
Еа ІП В; 
ПОТ КЕТИБ: 
zeros(4,1)]; 


%//end of file flatearthdyn.m 
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Fil у r ormdrag(A,alt,vel2 DSoost) 
$FORMDRAG Computes form drag for a missile with frontal 


% area A іп a standard atmosphere 

% FORMDRAG (A, ALT, VEL2 , BOOST) 

% uses MACHVALT, CDVMACH, RHOVALT 

% see also 

$ Copyright 1999-2000 by Triple B Enterprises 


фокс 
kkkkkkk kk 

%// File: formdrag.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

%// Date: 9 May 00 

%// Description: Computes form drag for a missile with frontal 
% area A in a standard atmosphere 

%// Inputs: area, altitude, V^2, boost on/off 

$// Outputs: parasitic drag force 

$// Process: 

$// Assumptions: 

$// Warnings: 


Sy fe RARER REAR EE EERE RRR KKK RRR RK RK RR RRR A a КИР У а ааа 
ххххххххх 


%// Order of elements 


%// -Define globals 

%// -Define constants 

$// -Define elements of input vector 
%// ~Functions 


Бр Кккк кх кх ххх ххх Ххх к ORE A KE AR ТАЖ СА E а лла а а ЗЕ 
* * í * * KKK 


В УШИ аеғтпе о1ораї15 “У 
в". спе солпесапса шаа 
р ЕЕ мессез 5500 
оо ва не уата ео“ 


rÉrhoernovalctialtor 
machswwel2)^(1/2)/machvalt(alt):; 


% жута за: «тв во =» functions хжхххххххххкхх 
ка “= (mac 100) 


тасла.83: 
епа 
Cdzcdvmacnimach,Doost):; 


во мена РСА А 


%//end of file formdrag.m 
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foncer A gravity (u) 

SGRAVITY Computes simple gravity model for 6DOF model 
% see also 

% Copyright 1999-2000 by Triple B Enterprises 


%5//ХХХХХХХХХХХХХХЕХХХХХХХХХХЖХХХХЕХХХХЖЖХХХХХХЖЖЖХЕХЖЖХХХХХХХХХХХЕХ ЖЖЖ ke e 
хжхххххххх 

$// File: gravity.m 

$// Name:  LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler:  MatLab v5.3 


%// Date: 19 бер 00 
$// Description: computes spherical earth gravity for TBM 


$// target dynamics 
%// Imputs: target state vector 
%// Outputs: gravity vector 

%// Process: 

%// Assumptions: 


%// Warnings: 
$ / / * * * * * * * * k < * ke e eee e ke ce k k k kk kk kk K KK k k kk Kk'k k kk k kk k kk kk kk k kk kk k kk kk kus 


жхххжххххх 


%// Order of elements 


% / / -Define globals 

$// -Define constants 

%// -Define elements of input vector 
%// -Functions 


% / / * he ce he eee koe oe he he e e eoe ke e he hehe e e ke ck ke e de ke oe ee ce e e e ee e ke CK Se ke ok ee ce e cock kk ck Sk ck ck ck k Kk 
ххххххххх 


= На“ које пе globals хз 


$ ****** define constants ****** 
СМ Е-3.9860014е14; % G*mass of earth 


Ф АИ ОЕ Vector “я 


ааа ге макгларјег ии 
Imaq=sqrt (u (bi Cru а а 62): 


% ххжхххххххххх functions хххххххххххх 


y= (GM E/mag шча 


%//end of file gravity.m 


function ysgravity2(u) 

%GRAVITY2 Computes simple gravity model for 6DOF model 
% see also 

Ф Copyright 1999-2000 by Triple B Enterprises 


ӘЛИ Ак жа сак O ккк К АКК КАК Кы АК cock cse ok e ok oe oko cec o X X A X AREER EEK REE KK Xe 
ххххххххх 

$// File: gravity2.m 

$// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$// Date: 19 бер 00 

%// Description: computes spherical earth gravity model 
Ф// for interceptor 

%// Inputs: missile state vector 

$// Outputs: gravity vector 

%// Process: 

%// Assumptions: 

%// Warnings: 


А о 
ххжхххххх 


%// Order of elements 


$// -Define globals 

$// -Define constants 

%// -Define elements of input vector 
%// -Functions 


кк OX КОМ ИЛИЈА КЛУ кх ккк ккк ххх КК АИ А coco Xe RK EUG 
* * x< xx x x kx 


фк Stine globals Sm = 


ф “ЖЕНА сене Constants “< 
СМ Е-3.9860014е14; % G*mass of earth 


Ф XXRAAK Serine Imput vector Хе 


5 ххх” initialize varlaplesS ““-- 
mag sar (ü(1)%2:0u(2)20 40 2): 


% хххххххххххх functions хххххххххххх 


yY- OMBRE mag 97 Титан T 


%//end of file gravity2.m 
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% / / * * k kk kk k k k kk k kk kk kk k k K K kk k Хх» ххх» +» RK kk kk KK Kk KK KK KK KK RK KK 


хххххжххх 


%// File: roe biol ER IH 

$// Маме: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler: MatLab v5.3 

%// Date: 7 Aug 00 

%// Description: Automatically computes a kinematic boundary using 
6DOF 


$// simulator with tracking Ете 
% / / -Streamlined search loops 

$// -Status indicator 

%// -Saves most recent data to disk 
% / / -Derived from KBOUTER 


%// Imputs: none 

%// Outputs: one figure of kinematic boundary 

%// Process: streamlined brute force search algorithm 
%$// Assumptions: попе 


$// Warnings: попе 
SJ f CK e echec eee ee Ke ee CK e e CU eS C Sec e OI COCOS Ke e ko eoe oec E 


ххххххххх 


%// Order of elements 


%// -Define globals 

Ф// -Define constants 

% / / -Define elements of input vector 
$// -Initialize variables 

$// -Functions 


qp] f kk kx kc kk Sk Ok КА КАКАН 


koc ck o cock ko Ж 


фо SS cetine globals "S559 


ф хе“ Sdetine Constants “ее 
thesisinit 

% set min engagement range (10000 m default) 
minrng=10000; 

$ set heading increment 

degstep=5; 


$’ -de ine тпрпкЕ Vector 555 
G tte ee Qe a Ze variables TAa 


Mast = Minnes Е 

load current 

% set target altitude 

СОСЕД EU $ default co-altitude 

$ set target turn rate. default=0 degrees/sec 
target _ turn- 0; 

% set target speed 

раштасл=н % user sets Mach # 
tgtspd=tgtmach*machvalt(tgtalt);% machine computes speed 
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% хххххххххххх functions хххххххххххх 

% start in tall chase step to head on by <degstep> 
% degree increments 

for heading=0:degstep:180 


cic 
plotcount=1; runplot=[]; rangemax=0; 
heading % show heading counter 


tgthdg-heading*pi/180; 


% compute target speed components 
xspd-tgtspd*cos(tgthdg); 
yspd=tgtspd*sin(tgthdg) ; 


% first range loop step by 10 km 

Сот  Egtrngsminrng: 1060060 1 0 OD 
алер [5 та БРЕ (еј Белој лан s As) 
% set initial target state 
саха Едва spo O spo аса ЕЕ; 
% initialize filter 


ALAST=(EGtinat (1) ;totinit (2c ini es) egenic Ae O CEES) л 
Sco) ОИ 


% call simulation 
Sim(’thesislfilt’) 


% analyze data from current run 

range-sqrt((MissileOut(:,1)-TgtOut(:,1)).^2*... 
(MisslleOut(:,2)-TgtOut(:,3)).“”2+... 
(Mise):leOnt(i:,3)-TgbOUt(ic 5)).^2): 

$ save run data 

тито те (Ро ве ои а] = | пад аме crj 

$ score run 

if (min(range)»5) 
break 

end 

plotcount=plotcount+1; 

end 


idx=find(runplot 2!” 

Lt (ise) 
rangemax=runplot (max (idx) ,2); 

end 

ва ја = |]: 

plotcount=1; 
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% 1 Кш step size for max range. Streamlining code 
tgtrng=rangemax+1000; 

di pPI FAA umumPscr(tgerng),' maxik ^ 21] 

$ set initial target state 
LogtgnmE-[|tgtrng;xsbd 0 wespd total: 0]. 

$ initialize filter 


xXbAST=[Egtinit (1)”7? CgCinit (2)%0; EgGE1nmit (3) САСА ОКИС Ене Е ЕЕ Я 
ЕКОО: 


€ call simulation 
sim(’thesislifilt’) 


$ analyze data from current run 
rangessqrt((MissileOut(:,l)-TgtOut(:,1)).^2-*... 
(MissileOut(:,2)-TqtOnt(05 оне ве 
(MissileQut(:,3)-TgtOUE| 255220 
1f (min(range)<=5) 
rangemax=tgtrng; 
tgtrng=tgtrng+4000; 
disp И Тоо Ест (саспа атк" "У 
% set initial target state 
ЕО ЕЕ ЗА = ЕО xsSpada U; yspd;tgtalt 0]; 
Ф initialize filter 


XLASTs[tgtinit(1);tgtznzt(2);0;tgtindtt(5);tgtintte(4); 0G ng 5) eo 
intense 


% call simulation 
sim(’thesislfilt’) 


$ analyze data from current run 
ranye=sdgxrt ((MiSSiTeOnE (: | у-тавоне о 
(MissileOut(:,2)-TgtOut(:,3)).^2-*... 
Мага а veOut (2473) -TotOut(:¢5)). 2) 
if (min(range) <=5) 
rangemax-tgtrng; 
end 
end 
$ main search loop ikm step size 
for tgtrng=rangemax+1000:1000: (rangemax+4000) 


disp( i *** T? mumZstr(totrna),. па ИЯ 
% set initial target state 
LGEIinicSe[torrug:;xSpa; 0 YSP; totalt Ol; 

% initialize filter 


XLASTs[tgtintrE(l1);tgEIDIt(2);0.tot dimid t3) totini л О Cge nre а 
ir (06907 


$ call simulation 
sim chesisliilt ) 
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$ analyze data from current run 

rancdcgessert((MissileOoup( - "uETOgtOUEt(:,1))-7^2* 5 
(MassileOut (:,2) =Tateur(: , 39mm 2+... 
(MissileOut(:5 3) TGpUNE(:. 599 ^2); 


$ save run data 
rEpumsLobiplotcount,')bom»nirsmge)tgtrng]s 
$ score run 
if (min(range)»5) 

break 
end 


piocEcount-plotcount-gzde 
end 


тах C nmd ЕШ НО dor 

Е 59, 
rangemax=runplot (max (idx) ,2); 

end 

runple s= [p] 

ploteount=1> 


$ 100 м step size for max range. Streamlined code. 
tgtrng=rangemax+100; 
а ери ј =>“ =" num2stet eqrerme x “песен ШТІ 
$ set initial target state 
pgEUDMESItqgtrngsxspdU vspe ЕС САЛЕ 
$ initialize filter 


XLASTe[tgtinit(1);tgtinit(2);0;tgtipuit (o Egtlimer (1) ЕСЕ ЕСЕН 
geo 


% call simulation 
sim(' thesis lr ilt’) 


% analyze data from current run 
гагде=5атЕ ( (Мі5511еопе(:, 1) =тасбОцЕ(: 1) ) 2+ 
(MyssqgfeQut(:;2) TgEOUNE:,3))5923. 05 
(MiSSileoOut(: 3 STg Om 5). 2); 
if (min(range) <=5) 
rangemax=tgtrng; 
tgtrngstgtrngs400:7 
arspi['***-” mum st ruscemnegee Нах100 1 
$ set initial target state 
сава с= [сасгла: =ра 0 Уу5ра саса ШЕ 
% initialize filter 


XDBAST-s[tgbiDniE(l);;tgberndit(2);0 Оты (л) Барпы a) Баета О EIS NNI 
Jr (Oo ире 


% call simulation 
sim(’thesislfilt’) 


$ analyze data from current run 

range-sqrt((MissileOut(:;l)-TgtOut DM 
(MissileOut(:,2)-mDOdgtOut(:,.3)) 9921.55 
(MissileOut(:,3) Tj OD C (0, 5)9 2); 

1Ё (min(range)<=5) 
rangemax=tgtrng; 

епа 

епа 


$ main search loop 100 m 
for tgtrng=rangemax+100:100: (rangemax+400) 


disp(['*** "> numZ2str(eotend) пас а ТІ 
% set initial target state 

EgeEiInit=(tgtrng; xspd-0>yspa-tetalt: 0]; 

% initialize filter 


XLASTs[tgtinit(l):tgtinit(2);0;tgtingt(3);tgtin2t( DI eot neno ЕС 
ТЕС ОЕ 


% call simulation 
sim Chesislrile 


$ analyze data from current run 
range=sqxrt((MissileOut(:,1)-TgtOut(:,1)).“°2+... 
(MissileOut(:,2)-TgtOut(:,3)). -2+... 
(Ма 55 11е0џ (:,3)-"тавов 5122 


$ save run data 
runplotiplotcount.:)-[maintrange)tgue- 
$ score run 
if (min(range)»5) 

break 
end 


plotcount=plotcount+1; 
end 


Ва = (типо (з, ==); 
пе (тахо) 
rangemax-runplot (max(idx),2); 
end 
runplot=[]; 
plotcount=1; 
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% 10 m step size for max range. Streamlined code 
tgtrng=rangemax+10; 

sea те век собе и Inacl)***']3 

% set initial target state 

еп Ее= ГЕЈЕ; хера ОгуграредјЕа12; 0); 

% initialize filter 


ОТ (КЕ о РЕЈА РОВ а (СА О ере Пана ере CgEInit (a ОЕЕО ЕЕЕ СЬ у Ея 
ЕО 


€ call simulation 
Simi tChesisifilt’) 


% analyze data from current кип 
range=sqrt ( (MissileOut(:,1)-TgtOut(:,1)).“ ”2+... 
(MissileOut(:,2)-TgtOuE(:;3)] 7 ^29 
(MrissleQut(:;3)-IgtOUE O Soe. 2). 
ТЕ (пап(веосе)<=5) 
rangemax=tgtrng; 
tgtrng=tgtrng+40; 
Слер (ү t a "што ври Гаспа азс ВО) 
% set initial target state 
EGtinue—(EGQErNG-xspa: 0 vend Eqrale. 0]; 
% initialize filter 


ALAST=[tgtinit (1) ;tgqtinit (2); 0-tgtinict(s) БУЕ ЕИО ЕЕЕ P raps 
ЕООД 


$ са11 simulation 
и (ЕВЕ СВЕТЕ” 


% analyze data from current run 

rangessart((MissileQut(:; l)-TygtoOUt(: ¿1)) 22. 
(MissrleQut(:,2)-TgtOut' "3 ) o e 
(Miss)leOubt(:,3)--TIgEODDC(.,35)9898 2). 

if (min(range) <=5) 
rangemax=tgtrng; 

end 

end 


$ main search loop 10m. Note, we are now computing the 
+ rull outpüt vector for each run. 
for tgtrng=rangemax: 10: (rangemax+40) 


drsp([^*** "num2strdstrng 9c терсе) сағаны 
% set initial target state 
tgtinit=[tgtrng:;xspa;:0;yspa:Egtalkt i ols 

% initialize filter 


XLAST=[tgtinit(1);cgtinit (2); 0; tgtinit (3) tgtinit(4) 0 corre ее 
IEG OJ; 


% call simulation 
sim(’thesislfilt’ ) 


ШЕ? 


$ analyze data from current run 

range=sqrt((Missileeuc(: , 1) -Ttoteut@, Ye 
(М1 $ 511е04Е (52) Тов. 
(Маа5511ебц  (:„/3)—Ттавоне 5) 2); 
t=MissileOut(:,14); 

index=find(range==min(range) ); 

А ве ІШ” 


& Compute cost function J=20re(tI) 2FIinteg( uU 2 200 

% and missile divert 

u2=(omegaout(:,1).*2+omegaout(:,2).%2); 

integral=0; 

for il=2:index 
integral-integral-(t(ii)-t(ii-1))*u2(ii-1); 

end 

J=20*min (range) *2+integral/1000; 


% save run data [miss dist,cost,divert,time,max range] 
runplot (plotcount, :)=[min(vange) ,J,1ntegraln ip ыу сЕ 


ЈЕ (min(range)>5) 
break 
end 


plotcounrsplotcountet: 
end 


ара — ен а е runphioE ( -)e5)€ 
ТЕ (са) 

rangemax-runplot (max(idx),5); 
end 


if (isempty(idx) ) 

maxnit (neading+r1, :)=[0 0 050 0]; 
else 

Maxhat (headings? : )ë=grunplot (max (1qx) I; 
end 


runplot=[]; 
piotcounmEsl: 


£ save data to disk 
save current maxhit 


COC 
% note for some guidance laws, the down step here 
% must be 2 or more----------------- | 
minrng=10000* (floor (rangemax/10000)-1); 
Lf Aman eng 50007 

Ааа = ОООО 
епа 


епа 
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$ plot it for me baby 
rholzmaxhit(1:degstep:181,5); 
је = пе пола ома (кло1)]: 
theta=180:degstep:360; 
theta=pi/180*theta; 
theta=[theta,-1*fliplr(theta)]‘; 
figure(5) 

polar(theta,rhol) 


2//end of file KBOUTER.m 


141 


5//ХХХХХХХХХХХХХХХХХЕХХХХХХХХХХХХХХХХЕХХХХХХХХХХХЕХХХХХХХХХХХХХХХХХХА ХАЖ 


* * * ko ko ko ko * 


%// Flle: KBOUTER2.m 

%// Мапе: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

%$// Date: 7 Aug 00 

%// Description: Automatically computes a kinematic boundary using 
6DOF 


%// simulator. 

% / / -Streamlined search loops 

% / / -Status indicator 

% / / -Saves most recent data to disk 
$// -Derived from KBOUTER 


$// Inputs: none 

$// Outputs: one figure of kinematic boundary 

%// Process: streamlined brute force search algorithm 
%// Assumptions: попе 


%$// Warnings: попе 
gb f f o koe ecce ke oe e e e e e he e he e e eee ce e eee eec khe he e e e ke khe cce e e ke cce e e e e ke e e e e ke e he e ek kk ke ke 


хжхххххххх 


%// Order of elements 


$ / / -Define globals 

5// -Define constants 

% / / -Define elements of input vector 
% / / -Initialize variables 

%// -Functions 


$ / f[ Kk e e X e he ke ke hehe e e he he e e e e ee % * he e e ee ee he X e ee e e e e eec e ce ee ke ecce ce e dece e e e e e ke ke e e € 


ххххххххх 


Ен сора ****** 


ы жж СЕстШетсопскапееа 55552 
thesisinit 

% set target tuxrn rate, default=0 
target_turn=0; 

% set min engagement range (10000 m default) 
таа = 0000 

% set heading increment 

degstep=5; 


ФО МАУ п деғтле опе уессот Уя 
$ ****** initialize variables ****** 


maxhit=1]> minnie: 
load current 


$ set target altitude 
асе: титс Ы $ default co-altitude 


$ set target speed 
toEmach-4995* $ user sets Mach £f 


tgtspd-tgtmach*machvalt(tgtalt);$ machine computes speed 
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% хжххххххххххх стопе Kx Kx x xxx x x xxx 

% start in tail chase step to head on by <degstep> 
% degree increments 

for heading=0:degstep:180 


tie 
pPplotcountc:=l; runplot—(]- = rangemax—0; 
heading % show heading counter 


tgthdg=heading*pi/180; 


$ compute target speed components 
xspd-tgtspd*cos(tgthdg); 
yspd=tgtspd*sin(tgthdg) ; 


$ first range loop step by 10 km 

Loreccrug-manrpesd00009150/0900 
disp([^***t " müom2str(botsng p b) 
$ set initial target state 
tgtsnpst-[totrbng.;xspe уара Б акаш II 


% call simulation 
ӨШІН Thesisl') 


% analyze data from current run 

range-sgrt((MissileQutí(.:,l)-TgtOut(: 13)).^224. 2 
(Missileout (32) -TO our C О o m 
(Missile0ut(: 3) - TgEDUE(: DDA; 

% save run data 

уипрјов (рјовсошпеЕ, ·)= ен trance) EogtEDng]s 

$ score run 

її (шпхїп(тапсе)>5) 
break 

епа 

р ОСО РОВ и 1 


= 


ena 


а век = 2 pd (sunplor (721) =5 
рә 1554 

rangemax=runplot (max(idx) ,2); 
end 
BUTS ОЕ, 
Ре 6 
$ 1 km step size for max range. Streamlining code 
ctrng=rangemax+1000; 
зет те ве Се а бы со Та ааа ын 

set initzal target state 

EQLIMIL=|totrng: xspda; Ul yspa:- totale. 01. 


oe () rt 


я саг ьт артоп 
sim(’Thesisi’) 
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% analyze data from current run 
range=sqrt((MissileOut(:,1)—TgtOus (:,1)) ..22+ 


ЕЁ 


(MissileOut(: 2) PgbtOSUt(:.3)).^2* 95 
(MissileOont (l 3)=STgEO5n SB ^2); 

(min (range) <=5) 

rangemax=tgtrng; 

tgtrng=tgtrng+4000; 

атр са т тте ст Тесеігпа), лозе ле s. Qp 
% set initial target state 
tgtinit=[Cgtrng;xspa;0;yspa косак ш 


$% Call simulat Ion 
sim('Thesisl') 


$ analyze data from current run 


rangessqgrt((MissileOnt(:.1)s5TgtOmE UID 22 we 


(MissileQuE ( 2) Tg OUE T 2. 
(Метео, Е. 
ЈЕ (min(range)<=5) 
rangemax=tgtrng; 
end 


епа 


% main search loop lkm step size 
for tgtrng=rangemax+1000:1000: (rangemax+4000) 


adispul’*** " num2str(togtrng) пева ~*~" |p) 
$ set initial target state 
жерк шор је — [Ое ә ета е ә есәее Гура: еј сеш 0] 


$ call simulation 
sim('Thesi1isil:“) 


% analyze data from current run 


range=sqre ( (Missi leCut(: 1) -TaetOut(= 71). 2h E 


(Missi leout(s, 2) -TOtCOUET +, 3). 024. 5 
(MisslleQut(-3)-TgtOut(:,5)9892)0* 


$ save run data 
fZuunploti(plorcounct.:)-Imonirange) 'cgerug]; 
$ score run 
if (miní(range)»5) 

break 
end 


plotcount-plotcountii: 


end 


аа oti n s=): 


LE 


(ШОХ) 
rangemax-runplot (max(idx),2); 


end 
Tunpleu 1: 
plotcount-zl; 
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% 100 m step size for max range. Streamlined code. 
tgtrng=rangemax+100; 

Sp * = опата век (са клетва те 1007 +“ ]) 
% set initial target state 

tgt inme- лава разора саса: 01: 


% call simulation 
©лш(”ТҺе51581/) 


% analyze data from current run 
range=sqrt ((MiSS1ile9Out(: 1) )=TgítOnt ( GPP022+ К 
(Missa leOQutr(:,2)—TqQtOurl= 35) а ен... 
(Missa leout 2.3) —TqtoOut (7 5)))2.2) с 
if (min(range)<=5) 
rangēmax=tťtgtrng:; 
tgtrng=tgtrrig+400; 
Са БОНА АНА опат Str во Бете | mays Оро") 
% set initial target state 
Бе = ес ст хора „О ;-yopa-tota 130]: 


$ call simulation 
sim('Thesisl') 


$ analyze data from current run 

rangessqrt((MissileOut(:,1)-TgtOut(:,1)).^2-*... 
(MissileOut(:,2)-TgtOut(:;3))«.^24:.5 
(MissileOut(5M3»PmgtOut(:;5)).^2)5 

if (min(range)<=5) 
rangemax-tgtrng; 

end 

end 


$ main search loop 100 m 
for tgtrng=rangemax+100:100: (rangemax+400) 


disel src ЕЕ (tqQermig), па 9 
$ set initial target state 
LEgbtrsnst-begtrng;xspds0svepdetebalt.0l 


$ call simulation 
Ssim('Thesisl^') 


$ analyze data from current run 

range-sqgrt((MissileOut(:,1)-TgtOut(:,1)). 
(Missi beOQut (422) =TgtoOut (475) ) 52 5... 
(Missi leOut (2,3) -—TogtOuc ts, Sy). 2): 


в 
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% save run data 
runplot (plotcounc r: )=(min (тапде) ‚сост ЕКЕ 
$ score run 
if (min(range)»5) 
break 
end 


piotcount-plotceceubEsd 
end 


idx=find(runplot т а-ы 
т} 
rangemax-runplot (max(idx),2); 
end 
rump loc =| |; 
plotcountszl; 
% 10 m step size for max range.  Streamlined code 
са кпа=капдетах+10; 
р ое |) 
% set initial target state 
ros irnit=[toryngg: xspa; 0: yspdgd: tgral О 


% call simulation 
sim('Thesisl') 


$ analyze data from current run 
range-sqgrtí((MissileOut(:,1)-TgcOut( иии o, m 
(Missi PeQutl (=) 2) -—Tatoue(:,.3)) 22. 
(MassiPteOur(- 73) -Taccue (275) )- 72); 
if (min(range)==5) 
rangemax-tgtrng; 
ОС ва Suc 0 
disptl алати ле (ере maxi =e 
$ set initial target state 
Саал ОЕ ЕП хер О узра. воза 05 


> Cal) simulation 
стт (ТТ лес ЫК 


% analyze data from current run 

талаё=<сагт((М1ївєт еш (с, у=тасОоп к (:, гу у о ны 
(Missi tleout (2) 2 y= 1oeOurt 2-3, 5 29 X 
(Матевајеони 3) ТЕ 0) 2016 

її (min (range) 5) 
rangemax=tgtrng; 

end 

ena 


]46 


$ main search loop 10 m. Note, we are now computing the 
$ full output vector for each run. 
for tgtrng=rangemax:10: (rangemax+40) 


displ 44a. mum2str(egtrnd),” maxl0*** 7) 
$ set initial target state 
Бејт е= „раскид. хора; ог узра ава је: 51: 


% call simulation 
sim('Thesisl') 


$ analyze data from current run 

range=sqrt ( (MissileOut(:,1)-TgtOut(:,1)).^2+... 
(Missīileðut(:,2)-TgtOuüut(:,3))- “2+... 
(Мес тесте ТОО 392); 
t=Missile0ut(:, 14); 

index-find(range--min(range)); 

ip-t(index(1)); 


$ compute cost function J=20*e(tf£) *2+integ(u%2) /200 

$ and missile divert 

u2-(omegaout(:,1).^2«0megaout(:,2) .^2); 

integral=0; 

отн 
integral-integral-(tí(ii)-t(ii-1))*u2(ii-1); 

end 

J-20*min(range)^2-«integral/1000; 


$ save run data [miss dist,cost,divert,time,max range] 
ranplot(plotcount,:)s[min(range) J5.a3utedgrsl ои Боен: 


if (mim (range)>5) 
break 
end 


püuocbtcount-plotcounttl: 
end 


Tax- -find про, ЕО): 
JE rds) 

rangemax-runplot (max(idx),5); 
end 


if (isempty(idx)) 


maxi; (peading+1 :)=[6 0.0 0 0]; 
else 

maxhit (heading-1,:)-2runplot(max(idx),:); 
end 
runplot-s[]: 


photcoss-l: 
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% save data to disk 
save current maxhit 
COC 
% note for some guidance laws, the down step here 
$ must be 2 or more----------------- | 
minrng=10000* (floor ( rangemax, ОООО = 
Lf (општој = 5000) 
minrng=5000; 
end 


end 


% РОК О ШЕ for me baby 
rhoismaxhut(l:degstep: T9075» 
тело, ЕТО (1001 
theta=180:degstep:360; 
theta=pi/180*theta; 
theta=[theta,-1*fliplr(theta)]’; 
figure(5) 

polar(theta,rhol) 


2£//end o£ file KBOUTER.m 
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function y=machvalt (alt) 
ФМАСНМАШТ Computes linear approximation for a given 


% altitude in meters/sec based on standard ICAO 
% atmosphere 

% MACHVALT(ALT) 

% see also CDVMACH 

% Copyright 1999-2000 by Triple B Enterprises 


E так ьо оо ааа 
kx x k k k k k X 

%// File: machvalt.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

$// Date: 8 Jun 00 

%// Description: computes linear approximation to Mach 1 for 
% / / standard ICAO atmosphere 

%// Inputs: altitude 

%// Outputs: Mach 1 

%// Process: 

%// Assumptions: 


%// Warnings: 
БАКА АКК КК k k k k X x < Kk ok oko k ok kock KOKOK OK к OX UXOR 


*ockockock k k k k xk 


2// Order of elements 


% / / -Define globals 

% / / -Define constants 

%// -Define elements of input vector 
$// -Functions 


B/ [RRR KKK RRR KKK KKK KK ERK KKK EKER KEK KKK KKK KEKKEKKEKEKEKEKKKKKKKEKEKKEKKKKKKK KEK 
ххххххххх 


БК мене ексера ке je **** <* 


р пе лсопе алке 75“ 
Machi=[-.0041 340.3]; 

Мнеп2-205--7- 

Меспз-| 09067 2871 


m *****x derine input vectors serm 


к от рлар те муа араз e > 
а р=аре (а1 5): % account for NED coords 


% хххххххххххх functions хххххххххххх 
ТЕ (асер) 
y=polyval(Machi,alt) ; 
else 
if (alt>20000) 
y=polyval(Mach3,alt); 
else 
y=Mach2; 
end 
end 


2//end of file machvalt.m 
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5//ХХЕХХХХЕХХХХЖХХХХХХХХХХХХХХХХХХХХХХЕХЕХХЕХХЕХЕХХХХХХХХХХХХХХХЕХХЕХЕХЖЕХЖЕезлЖ 
* kx k Ж 

$// File: missiledata.m 

%// Мапе:  LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 17 May 00 

%// Description: missile data for AMRAAM 

$// Inputs: попе 

%// Outputs: various 

%// Process: 

%// Assumptions: 


%// Warnings: 
Ф-//ХХЕКХК ХЕ КК Кик KK E 


ал жж жж Ж 


%// Order of elements 


% / / -Define globals 

%// -Define constants 

%// -Define elements of input vector 
$// -Functions 


фу A AEAEE E KEKE KR KER A ER A eee 
* k * * k x Kk x x 
% Establlshes missile dimensions for use in computing 


% Aerodynamic forces and moments 
% Except where noted, all dimensions in MKS system 


$ Missile Name: PSEUDO AMRAAM 

mut xc ш па тест о Ба еге зк 
global m d L XCG XCPN XCPW XCPB XHL 
global ST SW SPLAN SREF 


$ ****** define constants ****** 


%——— missile body dimensions ----------------------- 
ПЕ 56667 $ mass, may be time varying 

а= 1/76; % diameter 

Поло % length 

Хеб 82388; $ initial c.g., may be time varying 
LN-.6769; $ length of nose cone 

%------------ tailplane dimensions -------------------------- 
XHL=3.454; $ hinge line arm 

СЕТ=.4061; ж Бал: root слота 

C p= ОРУ. stail tip cCrhorad 

= O06 O; % tail extension 

D 22507 $ tail height 

%--------------- wing dimensions --------------------------- 
XW-z1.134; £ wing to radome tangency point 
CRW=.3554; $ wing root chord 

CASU iwing tip chord 

wW rO; % wing extension 

НИ Уа; % wing height 
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ТО КЕКЕ» ЯЕ re input vector 


ф кик Initialize variables 


*ockockockock ok 


kak kkk 


$ хххххххххххх functions * * * * ck ockock xx x 


ната centers of presse 


Е РИ=ЕСУ ЕМ 
XCPW=LN+XW+. 7*CRW-.2*CTW; 
AN=.67*LN*d; 
Ав=( а: 


XCPB=(.67*AN*LN+AB* (LN+.5* (L-LN)))/... 


(AN+AB) ; 


%-------------- area computation 


SW= .5*HW* (CTW+CRW) +CRW*WXT; 
Sis. >*HT*(CTT+CRT) +CRI АТ? 


% 
% 


wing area 
tail area 


> > о eee — жамы» 


wing CP 
plan area of nose 
plan area of body 
body CP 


| ақы» eee eee eee олы» бы» ақы» аңы» «шы» D = 


SPLAN= (L-LN) *d+.67*L*d; % body and nose plan area 
SREF=p1*d°*2/4; $ missile cross section 
%------------- compute the inertial matrix ------------------ 
т=а/2; 

dx-m*r^2/2: 


Чу - те рогу врати а Ајент 05/25 x GG) 527 


Jz=Jvy; 


%//епа of file missiledata.m 
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Qs f [ke ke hee he e ke e e ke he ke e ke ke e ke e hee he ehe ce ke ke e ke ce he ce e ce ke e ke ce ee ke ee ke e e he ce ke ce ck KR KK KK kK KR KK 
хжхххххххх 


%// File: missiledata2.m 

%// Name: LCDR Robert D. Broadston 
%// MSEE/EE Thesis 

5) Op rating Environment: 
%// Compiler: MatLab v5.3 
%// Date: 12 Арк 00 

%// Description: missile data for Jerger missile from Zarchan 
%// Inputs: none 

47 / mOUEDULS: various 

%// Process: 

%// Assumptions: 


%// Warnings: 
Z// k kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


Windows NT 4.0 Service Pack 5 


kk kk kkk kxk 


2// Order of elements 


$// -Define globals 

$// -Define constants 

$// -Define elements of input vector 
$// -Functions 


Z// kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
ххххххххх 


% Establishes missile dimensions for use in computing 
% Aerodynamic forces and moments 
% Except where noted, all dimensions in MKS system 


$ Missile Name:  JERGER 


$ хжххххх хжхххх 


define globals 
global m d L XCG XCPN XCPW XCPB XHL 
global ST SW SPLAN SREF 


Фф ****** define constants  *****x 


%------------ missile body dimensions ----------------------- 
ПАЗ А $ mass, may be time varying 

= БР $ diameter 

Eon elo. % length 

ХСС=3.048; $ initial c.g., may be time varying 
= оса" $ length of nose cone 

$———— tailplane dimensions -------------------------- 
ХНЬ=5.9436; % hinge line arm 

СЕТ=.6096; % tail тосо соста 

CTT=. 0; $ tail Ciprchere 

eee: % tail extension 

HT=.6096; $ tail height 

%--------------- wing dimensions --------------------------- 
21080 $ wing to radome tangency point 
CRW-1.8288; t wing root епота 

CTWz0; % wing tip cnord 

WAT=0- % wing extension 

HW=.6096; % wing height 


Ба Cle ome, INnpmee Vector. ****** 


5 кке Anni tialize Variables ****** 


9; * < *& ck kK k Kk K k + xk LUNCtions хххххххххххх 


centers of pressure 


cee mew = = = = = 


тәме «жы» «жие «жж» «жы» «шы» «ше» «шы» ---- = 


XEePN= СУШЕ $ nose CP 
XC PW=LN+XW+.7*CRW-.2*CTW; $ wing CP 
AN=.67*LN*d; $ plan area of nose 
AB- (L-LN) *d; $ plan area of body 
XCPB=(.67*AN*LN+AB* (LN+.5* (L-LN)))/... % body СР 

(AN+AB) ; 
%-------------- area computation ---------------------------- 
SW= .5*HW* (CTW+CRW) +CRW*WxT; $ wing area 
S= БРАН (CTETCRT)+GRE“ TXT; % tall area 
SPLAN- (L-LN) *d«.67*L*d; $ body and nose plan area 
SREF=pi*d*2/4; % missile cross section 


%----------- compute the inertial matrix 
r=d/2; 

Пе 0200: 

ОУ УОТА 

Jz=Jy; 


%//end of file missiledata2.m 


"ыш ғ» қы қыш «жы» «лы» «шы» «шы» «шы» «қыш «шы» «қыш «шы» «қыш — 


S//**k***x* 
жхххххххх 
$// File: 
%// Name: 
%// 
%// 
% / / 
%// 
%// 
$// 
$// 
$// 
%// 
%// 
S//*****x 
kkkkkkkk k 
%// Order 
%// 

%// 

%// 

% / / 

ки 
хжххххххх 
% Establi 
% 


Орега 
Compl 
Date: 
Descr 
Input 
Qutpu 
Proce 


$ Except where noted, 


$ Missile 


хжжххх 


$ 


хххххххххххххххххххххххххххххххххххххххххххххххххххххххххххххх 


missiledata3.m 
LCDR Robert D. Broadston 


MSEE/EE Thesis 


ting Environment: Windows NT 4.0 Service Pack 5 
ler: MatLab v5.3 
25 Aug 00 
iption: computes missile data for SM-2 MR 
S: попе 
tS: various 


SS: 


Assumptions: 
Warnings: 


KKK KKK KEK KKK KKK KKK KKK KKK KKK EK KKK KKK KKK KK xx x xx x x x kk kx x k kx 


of elements 

-Define globals 

-Define constants 

-Define elements of input vector 


-Functions 
ххжххххххххххххххххххххххххххххххххххкхххххххххххххххххххххххххххх 


shes missile dimensions for use in computing 


Aerodynamic forces and moments 


all dimensions in MKS system 


Name: STANDARD RIM-67 MR 


*ckock ko k 


define globals 


global m d L XCG XCPN XCPW XCPB XHL 


global ST 


*ockokock k X 


xOGIO 2005 
EN=. 728 
= ees oe 
СЕЕ 3506; 
те 127; 
О. 
НТЕ, 
Ше! 12- 
FI 2 314 
Cu =] 93; 
WXTzO: 

Ru 42; 


SW SPLAN SREF 


define constants  ****** 


missile body dimensions 


$ mass, may be time varying 
$ diameter 
$ length 
; $ initial c.g., may be time varying 


$ length of nose cone 
tailplane dimensions 
% hinge line arm 
$ tail root chord 
t ral tip chord 
% tall extension 
% tall height 
wlng dimensions 
% wing to radome tangency point 
% wing root chord 
+ wing tip chord 
% wing extension 
% wing height 


— — = = — — — чы — — — «жы» «лы» «шы» «шы» «шы» «шы» «шы» «шы» «шы» «шы» «шы» «шы» «шы» «шы» олы» «шы» [ т== 


, 


= — --- --- --- --- --- — ешю «ши» «ли» «лы» ели» «ши» «шы» -... .... «ши» -.. «низ «жы» «шы»  — шш шш =—— —=— — — 


, 
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ы сегіпе input veetor  *****x* 


% ****** jnitlallze variables  ****** 


% хххххххххххх tunectrone ххх 
%-------------- centers of pressure ------------------------- 
XCEN=. 677 LN; % nose CP 
XCPW=LN+XW+.7*CRW-.2*CTW; % wing CP 
AN=.67*LN#*d; % plan area of nose 
ABz (L-LN) *d; $ plan area of body 
XCPB=(.67*AN*LN+AB* (LN+.5*(L-LN)))/... % body CP 

(AN+AB) ; 
%-------------- area computation ---------------------------- 
SwW=.5*HW* (CTW+CRW) +CRW*WXT; $ wing area 
5Т=.5*НТ* (СТТ+СЕТ) +САТ*ТХТ; % tail area 
SPLAN= (L-LN) *d+.67*L*d; $ body and nose plan area 


оо 


ОЗ ее 1: missile cross section 


%----------- compute the inertial matrix -------------------- 
ey 2; 

dorem*r 2642 

Syem* (bc 2/ i 2+n-2/ 4d): 

JZ=Jy; 


¢//end of file missiledata3.m 
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S5 f f kk ok e ke ke oe K X k Y X e ke X X ke eee e XK X X хх хх хх хх e hee e ke ee ee ЧЕККЕ ce kk he ck kk e ke ke ke 
KRKKKKKKKK 


%// File: missiledata4.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 19 Sep 00 

%// Description: computes missile data for SM-2 ER 
%// Inputs: попе 

%// Outputs: none 

%// Process: 

%// Assumptions: 

%// Warnings: 


$] у «Жк ЖЖ KX Xe 
oco cock cA ck kkk 


%// Order of elements 


%// -Define globals 

%// -Define constants 

%// -Define elements of input vector 
$// -Functions 


фу у «ЖЖ ККИ КК КК He 
ххххххххх 
% Establishes missile dimensions for use in computing 


% Aerodynamic forces and moments 
% Except where noted, all dimensions in MKS system 


% Missile Name: STANDARD WITH BOOSTER 
фо derine glopals SFr 
global m d L XCG XCPN XCPW XCPB XHL 
global ST SW SPLAN SREF 


% ****** define constants ****** 


%------------ missile body dimensions ----------------------- 
mn2d 5e $ mass, may be time varying 

ее, % diameter 
1=7.976; % length 
XCG=3.967; % Initial c.g., may be time varying 
ИН. 729: $ length of nose cone 

%———— tallplane dimensions -------------------------- 
Сб % hinge line arm 

СВТ. % tail root chord 

(Te 3- Ф Бат гір chord 
TXT: 07 % tail extension 

ES OE $ tail height 

%--------------- wing dimensions --------------------------- 
ке О; % wing to radome tangency point 
CRW-2.7; % wing root Chord 

CTW=2.5; Ss wing tip chord 
WAT=0; % wing extension 

HW=.151; $ wing height 
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с input vector ****** 


о хх _толеша те маглартев == === 


% хххххххххххх functions хххххххххххх 
%-------------- centers of pressure ------------------------- 
XS tme $ nose CP 
XCPW=LN+XW+.7*CRW-.2*CTW; $ wing CP 
РО LN a: % plan area of nose 
AB- (L-LN) *d; $ plan area of body 
XCPB2(.67*AN*LN-*AB* (LN* .5* (L-LN) ) )/... $ body CP 

(AN+AB) ; 
%-------------- area computation ---------------------------- 
SW=.5*HW* (CTW+CRW) +CRW*WXT; $ wing area 
ST=.5*HT* (CTT+CRT)+CRT*TXT; % tail area 
SPLAN= (L-LN) *d+.67*L*d; $ body and nose plan area 
SREF-pi*d^2/4; $ missile cross section 


%----------- compute the inertial matrix -------------------- 
r=d/2; 

а 

ОЕ ЕТО ОХА) 

Јг=Ју; 


%//end of file missiledata4.m 
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function y=modelswitch(u) 
$MODELSWITCH Switches missile models to simulate 


% staging for TBM interceptor demo 

% MODELSWITCH (T) 

% see also 

% Copyright 1999-2000 bv Triple B Enterprises 


% / / * * kk k k K T e e hehe ee ecce cese ke ke cese e e ke ke cce ke kk ke e c e ke ke ke e ke ck ke ke e ck ck ck ck ck ck ck ck ko ke ke Sk ck ko ke kk 


*ckckocko ck k ck k 


$// File: modelswitch.m 

$// Name:  LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

%// расе: 19 бер 2000 

%// Description: Switches between interceptor with booster 
57 and without booster for TBM demo 

%// Inputs: simulation time 

$// Outputs: попе 

$// Process: none 

$// Assumptions: 

$// Warnings: 


8b [RC koe ck e kk he e e kk ke ke Sk e ke e ck ee ce ke ke ck ck ck ck ck ke ko ce ck Cocks ck ck ck ck ck ck ck ck ck ck Ck ck ck ck Ck ck ck ck ck ck ck ck kk Ck ko ko 


Kk x x wx x x x x 


2// Order of elements 


$// -Define globals 

% / / -Define constants 

2/7 -Define elements of input vector 
% / / -Functions 


SS f f[ Rok o eee ee e ke Se Sk ke he Sk e se ee eee e khe kk ehe cce K e ce de ke ke se ck he ck e ck ck ke ck ck kk ock hock KKK KKKA 


x xxx А ХК ЖЕ 


а 
% а Ее Е > 
во лек“  дегтплпе Input vector 74505 
ФЕ ҰЗ initialize variables:  ****** 
+ “wee ere eee eee F Емас ЕТ оне хххххххххххх 
ЕО) 
missilecata3; 
else 
missileqaca4; 


end 


%// end of file modelswitch.m 
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жж жак КЖ ЖАКЕ К Ж жох хх Хх Хх хк Ккк кк ху хакк ккк К кк 
Хххххж ххх 


$// File: noisestudy.m 

$// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%$// Date: 19 Sep 00 

%// Description: performs noise study using modified AAM model 
%// thesisnoise.mdl 

$// Inputs: 

$// Qutputs: 

$// Process: 

$// Assumptions: 


$// Warnings: 
%// ХХХ КА КК KEE KEKE KRKEEKEERER EKER EERE KERR ER KK EK KKK KKK K 


ххххххххх 


%// Order of elements 


%// -Define globals 

$// -Define constants 

%// -Define elements of input vector 
5 -Functions 


b Рк А АА ААҚ а ААА АЖЖ ЖЖ АА ЖАК Ж ЖКС Ж Ж ЖА АСА ЖАСАСА 2 Z А а 
kek KKK KKK 


по А 5 
io efine Cons ants ~ <= 
па" Уаестежт риш ес Еос t tr x tx 


poc nmgtialize variables ****** 
$ initialize simulation 

thesisinit 

$ initialize variables 

holdrange-[]; 

holdpos=[]; 

% initialize target 

(ера == јеге о (122190.6000. 25): 

ЈЕ 
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% хжхххххххххх functions хххххххххххх 


% 100 realizations 
for numloops=1:100 
disp (numloops) 
sim('thesisnoise') 
$ analyze data from current run 
rangessqrt((MissileOut(:); b-TgtOUt(:,1)) 291 == 
(MissileOut (:,2)-TgtOut(:,3)).^2+.. 
(Misslle9out(:;,3) SR; Ош (: ,5)). 22), 
disp(min(range) ) 
holdrange (numloops) =min (range) ; 
1dx=find(range==min (range) ); 
hnoldpos (numloops,..)=Missileouciidx, 1:3) -TatOuc (1a) 2 
end 
missdistance=mean (holdrange) 
sigmadistance-std(holdrange) 


figure(5) 
plot (boiapos(: |) nolapaos(: ,2) ,Holgeos(: ји а а) 


$//end of file noisestudy.m 
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function v propnhayv3D(u) 
SPROPNAV3D Computes exact proportional navigation 


% in three dimensions for full aero model 
% see also EXACTPROPNAV2 
% Copyright 1999-2000 by Triple B Enterprises 


БТА та АД АЕ АЕ ЈАВА СА ДЕЛ ЈЕ АЕ ОВЕ ЛЕ REOR UK OK X ERU А ЗА т атта 
ххххххххх 

T Eile: Dropnav3d.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 12 Sep 00 

%// Description: 3D prop nav law for full aero model 

%// Inputs: seeker output 

%// Outputs: command accelerations 

%// Process: 

%// Assumptions: 


%// Warnings: 
EA ххх ххх ххх к хх А А АА АККА КАК КА OE UK US OE OK xc esc c eoo ic coke SKK cox oe ceo oko oc o oe EUR ooo: 


oko ck ok ck ok kk 


%// Order of elements 


%// -Define globals 

%// -Define constants 

%// -Define elements of input vector 
%// -Functions 


ng sex c RR ee 
хххжххххх 


ей 
global satflag 


ы 77575 define constants -***»** 
Nprime=5; 
Nprimez=5; 


о а прое <= --= 
thetadot=u(1); 
РЕНИ). 
loss) 
pnilos=n (a) 
Rsu(6); 

МС== (5); 
ћеад1лпа=ц (7); 
ла OU 
Vmdot=u (9); 
=). 
theta=u(11); 
psr; 
time=zu (13); 
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$ ******  initialize variables  ****** 


€ хжхххххххххх functions * * *k +) KKK Ke Kk Kk KK 
ny=Nprime*Vc* (thetadot) /cos(psi-los) ; 
nzsNprimez*Ve*tpbphidot)-9 75045». 


$% control force пісе 
if satflag 
if (abs(ny)»30*9.8045) 
ие 305922005; 
end 
ЈЕ (abs (nz)>30*9.8045) 
п2=519п(п2) *30* 9998045; 
епа 
епа 


Y= [nym |; 


$//end of file propnav3d.m 
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tunctdonew-probpnavptiu) 

%PROPNAVPT Computes exact proportional navigation 

% with dragforce inputs for point mass simulation 
see also EXACTPROPNAV2 

Copyright 1999-2000 by Triple B Enterprises 


oo ge 


БУЛК ххх хк куку ххх ххх ккк ккк X X ZX X X X X X X K K k ZR ОКО ОКА ОЉА ЖК СА 


хжхххххххх 


%// File: PROPNAVPT.m 

%// Мапе: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

$// Date: 24 May 2000 

%// Description: Proportional navigation guidance law for 6DOF 


% / / flight model. Computes applied forces for use 
% / / by induced drag model. Required to eliminate 
®/ / algebraic loops in the simulation 

t// ара [seeker data,IMU data,timer] 

%// Outputs: [command accelerations,applied forces] 


%// Process: proportional navigation law 
%// Assumptions: none 


%// Warnings: попе 
К ККК кк Кккк кк ERE LEER AERA EE EX ERR OE Ee Ee КАСА 


*oXkocko ko ko ok ok k k 


%// Order of elements 


%// -Define globals 

%// -Define constants 

% / / -Define elements of input vector 
%// -Functions 


аад 
*k* x *K e kx kx xk x 


2c detrne globajs 5 
global m satflag 


$ ****** define constants  ****** 
Nprime-5; 
Norimez=5; 


totes аегтпеті прое Vector ****** 
thetadotzu(1); 
phiGor u (2); 
Ше = o) 
pillos=u( (4); 
R=11(6) ; 

Ve “Sp (5); 
heading=u(7); 
Vm=u (8) ; 

Ма = (У је 
ла = о (100: 
theta=u(11); 
рѕі=0(12); 
Е1те=о (13); 
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% ******  initialize variables ****** 


% жжжххххххххх ппс Ес хжжххххххххх 

% classic PN guidance law 

ny=Nprime*Vc* (thetadot) /cos(psi-los) ; 

$ vertical acceleration must account for gravity 
nz-Nprimez*Vc*(phidot)/cos(theta-philos)-9.8045; 


= control force limiter 
if satflag 
if (abs (ny) >30*9. 38045} 
nveslgní(nv) 305903090257 
end 
if (abs(nz)>30*9.8045) 
ПЕ ВО) "30 S 6025. 
епа 
епа 


$ compute ABC forces applied 
Fx=0; 

Fy-ny*m; 

Е2-П2 “ІП; 


% output vector 
yv пара па FX EY EZ]; 


%//епа О file PROPNAVPT.m 
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tuncidconewespropmnavtlbmu) 
$PROPNAVTBM Computes exact proportional navigation 


% with dragforce inputs for point mass simulation 
% see also EXACTPROPNAV2 
% Copyright 1999-2000 by Triple B Enterprises 


% / / * * * *& * e ke e ehe k e k e k k k k e k k k k k kO kO K kO kO kO kO kO Kk KO kO e oe e e e KOK K K Ko KO KOK KOK RO K KO Ko KOR KO KO KOK KOK KO ЖК КЖК 
*x *< *K & < * kok * 


$// File: PROPNAVTBM.m 

$// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

$// Compiler: MatLab v5.3 

%// Date: 24 May 2000 

$// Description: Proportional navigation guidance law for 6DOF 


%// flight model. Computes applied forces for use 
%// by induced drag model. Required to eliminate 
$// algebraic loops in the simulation 

%// Inputs: [seeker data,IMU data, timer] 

$// Outputs: [command accelerations,applied forces] 


$// Process: proportional navigation law 
%// Assumptions: none 


$// Warnings: none 
ИХ ЖАК СА АА АТЖ ЖАЖА ЭЖЖ ЖА Ж ЭЖА ЖЖ Ж ЖАСА ҚАҚ Та АНА ane А 


KaKKK KKK KEK 


%// Order of elements 


%// -Define globals 

%// -Define constants 

%// -Define elements of input vector 
% / / -Functions 


Ж OR UKEOKUK KOKOK OK RUK ONE ORGOKOK OE OK OROREUK UN WORUR OK KO k kok k OK K K K K k kok k KOK kk K k'ik k Kk k SA kS eee 
ЖЖЖ ЖЖ ЖЖЖЖ 


na define globals ****** 
global m satflag 


$ ****** define constants ****** 
Nprime=5; 
Nprimez=5; 


Se etine imput vector ~****** 
thetadoc=u(1); 
Је = (2): 
los=u(3); 

phi los=1(4)- 
Е=и (6); 
Vc=-u(5) ; 
heading=u(7); 
Vm=u (8) ; 

еее to) 
posu qoe 
thera-uc br 
psi-u(12); 
time-zu(13); 
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Ф хи initialize уал зс c... 


% хххххххххххх functions хххххххххххх 


% classic PN guidance law 
ny=Nprime*Vc* (thetadot) ; 
$ vertical acceleration 
U2-NDIImez*Vvo*(phidobr)r 


$ control force limiter 
if satflag 
if (abs(ny)»30*9.8045) 
nvessrgnihy 3059 0456 
end 
if (absí(nz)»30*9.8045) 
nz=sign(nz)*%*30*9 6045; 
end 
end 


% compute ABC forces applied 
Fx=0; 

Fy=ny*m; 

Fz-nz*m; 


$ output vector 
yY nz AZ Fx Fy Fz]; 


%//end of file propnavtbm.m 
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function y=q2euler(u) 

S$Q2EULER Computes Euler angles from quaternions 

% see also ALPHABETA 

% Copyright 1999-2000 by Triple B Enterprises 


gs f f oko А АЕ А О О О КАКАО КАКО КОЈОЈ КАО АОЛ 


ok oko ck ck ko ok ok xk 


%// File: q2euler.m 

%// Мапе: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: .Matbabwv5.3 

%// Date: 6 Apr 00 

$// Description: computes euler angles from quaternions 
%// Inputs: quaternion 

$// Outputs: euler angles 

%// Process: Kuiper 

%// Assumptions: 

%// Warnings: 


Б/ ХХХ ККЖ ХХХ КАК ЖС ХЕ ЕЖ ЖЕЖ ЖЖК КИ т КЖ 
kak kkk KKK 


%// Order of elements 


%// -Define globals 

%// -Define constants 

% / / -Define elements of input vector 
% / / -Functions 


О Жл хек кк уух K X k X X k K X K X X KOK ххх K X X K ORO OK OK OK ORE жак кл атта аа а 
ххххххххх 


в ее с а15 ****** 

. derine соповатко <“ -<« 

Ё сарт пе прше vector w*****x 
а0-ч(1); 

gl=u(2); 

q2=u(3); 

q3=u(4); 


ВЕ и padeaalize variables , = 
% Ха Ах КЖ жж tunctrone ххжхххххххххх 
% convert quaternions to Euler angles 
mel—2*q0 2+2*ql-2-1; 

= ао“. 
mniss2*gitgs-2*g0*5402; 
о а и, 
03352*40092*2*0992—15 


psi=atan2 (ml2 mill); 

theta-asin(-ml3); 

$ о Е for singularity in pitch 

if (isreal(theta)) 
theta-theta; 

else 
thetassigni-mi5)*D1i72:; 

end 

phi=atan2 (m23,m33) ; 

y=([phi,theta,psi]J; 

$//end of file q2euler.m 
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ЕНЕНСЕШОП У =апас2р (67 
$OUAT2B Computes rotation matrix from quaternions 


% ОЏАТ2В (У) 
% see also QUATERNION, B2QUAT 
% Copyright 1999-2000 by Triple B Enterprises 


% / / * * * * * * e e ke k k k k k kk Ok k k kok K KOK K кх KO kO KO k Ko kR kO kO e KO ko kO kG K KO e KOK e se e e ck e ke kO kO R Ж 
ok ck ook oko ko ok ox 

$// File:  quat2b.m 

%// Мапе: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 8 Dec 99 

%// Description: computes rotation matrix from quaternions 
%// Inputs: quaternion 

$// Outputs: rotation matrix B 

$// Process:  Kuiper 

$// Assumptions: 


$// Warnings: 
Ф-//ХХЕХХЕХХЕХЕХХЕХХЕХХЕХЕХХХХХЕХХХХХХХЕЖХЕХЕХХЕТХХЕХЖХХХХЕХХХХХХХХХХХЖАСААВА 


хх > 


%// Order of elements 


$// -Define globals 

% / / -Define constants 

%// -Define elements of input vector 
5/7 -Functions 


А EEE REE RE AER ROARK RRR RAK ER A Oe eK MN Oe ee 
ххжхххххх 


о! с: 
толь > certinecconstdDUS ON не 


+ “== ЕСЕП vec Cor Er s s. 
а0=у (1); 

gqi-y (2); 
ату)” 
q3=y (4); 


$ ****** initialize variables  ****** 


% x kk КО * k * k * * functions хххххххххххх 


y=[q0^2+q1^2-q2^2-q3^2, 2*(а1*а2+а0*аз3з), 2* (ага аи 
2*(G1*G2=G0*G5)5 а0^2-а1^2+а2^2-а3^2, 2*(G2*G атта! 
092), 2 => (се "де ча), g0*2=-ql”" 2442" 2- G2 421 


%//end of file quat2b.m 
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function y=quaternion (phi,theta,psi) 


$SOUATERNION Computes quaternions from Euler angles 
% QUATERNION (PHI, THETA, PSI) 

% see also B2QUAT 

% Copyright 1999-2000 by Triple B Enterprises 


EA ko oe eoe oed ce e ce ce ke e e he e ke e ke ede ede e e ce ee echec ecco e echec ke oe ese ccce ck ehe e e ccce ce ck e ke e ke ke kk: 
* * * < * *< * kk 

%// File: quaternion.m 

%// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$// Date: 8 Aug 00 

$// Description: computes quaternion from euler angles 
%// Inputs: euler angles 

$// Outputs: quaternion 

$// Process:  Kuiper 

$// Assumptions: 

$// Warnings: 


Улу Кк кх клх клу ккк кх к хх к X X X k X X X K XX X K X K X X X K oko ook ook ин И А 
ххххххххх 


%// Order of elements 


%// -Define globals 

%// -Define constants 

$// -Define elements of input vector 
$// ~Functions 


Bf [= RRR RE K KIKIY RRR RR KR RK KK eR roc Kool ооо 
ххжхххххх 


аа ияи 

и «но decime constants. тя 

И сете input Vector —4****s 
клк Inalialize variables TERET 

с Евев Dy 172) *cCosSrtheta/2)*cosS(5sS17/2)5  .. 


sinlphi/2)*sin(theta/2)*siñ(psSi/2); 


aSr 05372) cos (Eheta/2)*cos( DsS172)= . >. 
SOS (Dhis2) *sin(thetay 2) *sini(psi72)- 


d2=CoS (Phi 2) sin (theta, 2) “cos (psi/2)442= 
Sin (omy 2) "cos (thetay 2) “sin (psi/ 2); 


q Cos e772) cos (theta/Z) *sin(psi/ 2)... 
Simipn./2)"sin(theta/2) “cos (psi/ 2); 


% * x * x * * кв ухкх огеле хжхххххххххх 


vesc epe oe rs 


$//end of file quaternion.m 
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Tunet ion y- rhovale(alt) 
SRHOVALT Computes atmospheric density vs altitude 


% for ICAO standard atmosphere 

% RHOVALT (ALT) 

% see also MACHVALT, CDVMACH 

% Copyright 1999-2000 by Triple B Enterprises 


фу ЖКХ кх ккк КККК КККК АК как жж кта ы 
жжххххххх 

%// File: rhovalt.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

%// Date: 8 Jun 00 

%// Description: computes atmospheric density from ICAO standard 
%// atmosphere. exponential model 

%// Inputs: altitude 

%// OWEDUES: rho 

%// Process: 

%// Assumptions: 


$// Warnings: 
Ф-//ХХХХЕХХХХХХХХЕАХ ХЖХЕХХХХХХХХЕХХХХЖКХЕХ ХХХ ck Sk kk e koe ck occ cR 


ххжххххххх 


%// Order of elements 


$// -Define globals 

%// -Define constants 

% / / -Define elements of input vector 
$// -Functions 


5//ХХХХХХХХХХЕЖХХЕХХХХХХХХХХХХХХХХХХХХЖЕХЖАХХХХХХХХАХХХХХХЖХАХХХХХХХХХХХЖАЖАЫ 
жхххххххх 


Е rs. rx; 
$ ****** define constants ****** 
в ее РЕ vector e a 


% ****** Initialize variables  ****** 
alt=abs(alt); % account for NED coord 


% * * * * *& * A k x Ax * * Оле топе *хх*ххХ% kx xx x Ж 


if alt>9144 

y=1 '75228763*exp (2>alr 5705. 6): 
else 

y=1.22557*exp(-alt/9144); 
end 


2//end o£ file rhovalt.m 
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БИНЕЕЛОП y=s1xdqordyn(u) 
$SIXDOFDYN Computes motion dymanics for six degrees 


$ of freedom 
% see also FLATEARTHDYN 
% Copyright 1999-2000 by Triple B Enterprises 


Акелла 
KKKKKKKKK 

$// File: sixdofdyn.m 

$// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 7 Арг 00 

%// Description: computes 6DOF dynamics in ECI coordinates 
%// Inputs: state vector 

$// Outputs: derivative of state vector 

$// Process: Stevens & Lewis 

%// Assumptions: 


$// Warnings: 
КУЛК а ААА Ж Ж ЖОЖ OK CK UK OK Ok KU ORO OK OK UK X Z XX obo Woo ook oe XO АСОВА А СО О 


ххххххххх 


%// Order of elements 


% / / —Define globals 

$ / / -Define constants 

%// —Define elements of input vector 
%// -Functions 


ff Re RAR RRR A NX RRR RRR XK KRONE Ko URN EUN Ko ee 
KKEKKKKKK 


о eer ine ајрора бака еже 


Ф ****** define constants ****** 


omega_x=7.292115e-5; $ earth rotation rate 

GM E-23.9860014e14; % G*mass of earth 

y E=6.375164e6; % radius of earth 

СЕ 982257; % ellipsoidal squash factor 


ae М ле ле Input vector  **"**** 
ыы): 

До = ИВ a(S); 66) J5 
omega b-[u(7);u(8) ;u(9)]; 
Pzu(7);  Qzu(8); R=u(9) 


f 


= (a GO) = eee ТАЗ) 4 
magq=norm(q,2); 
q=q/magq; 


х=[р;“_РБ:опеда_Б;а ); 
% inertial matrix 
J=[u(14) 0,00; 
Orga А O; 
0 0 o9 C P6)-]5 


% forces 
FOB= 01/17 ен а С 19091: 
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% torques 


PT Вв= и (20) о (21) а (22006 


$ the ever lovin’ force of gravity 
G p=[u(23);u(232)5 zea 


% and lest we forget, mass 


Шееле” 


% * * k * * x 


initialize variables 222: 
omega E-[omega x;0;0]; 


OMEGA, E-[0,0,0;0,0,-omega x;0,omega x,0]; 


$ compute rotation matrices 


B=quat2b(q); 


ОМЕСА_В=[0 -К 


К 
Б 


ОМЕСА а= [0 
-P 
=o 
-R 


% ххххххжхххж жж 


0 
p 


E 
0 
К 


-Q 


О; 
-Р; 

ШЕ 

O R; 
SINE 
0 -P; 
P 10 


functions kk Ok kk KK Kk KK 


% earth rotational velocity vector 


y-[OMEGA E, E zeros(3), Zeros (3745 
-B*OMEGA E^2, -(OMEGA, B-B*OMEGA E*B'), zeros(3), zeros (3,495 
zeros(3), zeros(3), -l*inv(J)*OMEGA B*J, zeros(3,4); 
zeros(4,3), zeros (4,3), zeros(4,3), -(1/2)*OMEGA q]; 


y-y*x; 


y=y+[zeros (3,1); 
Bto p+ F/m) HB; 


inv (J) TEB, 


zeros(4,1)]; 


$//end of file sixdofdyn.m 


фу * k k k * K k k k * k * k * * * k k k k k k k *k he e e e e k k k k k k k k k * Kk k k e ke ke К ИИС КК ККК ЖК ККК Ж К К 


ххххххххх 


%// 
SUN 
E 
> 
5 
%// 
%// 
БИ 
%// 
%// 
e7 
S 
uy 


File:  spielberg.m 

Name:  LCDR Robert D. Broadston 

MSEE/EE Thesis 

Operating Environment: Windows NT 4.0 Service Pack 5 

Compiler:  MatLab v5.3 

Date: 20 Sep 00 

Description: creates a movie of AAM simulation using 
thesism 

Inputs: none 

Outputs: none 

Process: 

Assumptions: 

Warnings: 


КК кк а 


ock * * * xx xx 


БИ 
%// 
А 
b 
БИ 


Order of elements 
-Define globals 
-Define constants 
-Define elements of input vector 
-Functions 


RL RRR ERR Ак Re KO ОТА РА А АДА ЕДИ ИО АЦЕ 


*ock ko kok ko kk 


Dux KW етен сера 5 екх 


$ ****** define constants ****** 


Б таене тпру есу 57 


$ ****** Jnitialize variables ****** 
thesisinit 


% *x* * k k * &K k k kx KKK fFuncelons KIRK KR RK KK KK 


figure(1) 


CUT 


Ба Бапа == <е0 2525 (40000, 5000,135); 
target_turn=12; 


for timestep=1:44 


tmax=.25* (timestep+0) 
sim('thesism') 
l=length(MissileOut) ; 

$ a little 3D action for the fans 
figure (1) 

С ЕГ 

view(-30,25) 

holdon 

axis([0 40000 0 8000 4500 755003) 


з 


% plot commands 
plercso (Missi1ileOut(:,;1) Miss1ileOQut (; 2)W=MISS ебе о јен 
Egtout(:,l); Таоци | -TgrtOUut(. 08 


plot3(MissileOut(1,1),MissileOut(1,2),-MissileOut(1,3),'*',... 
TgtOut (1,1), Tot@mimil; 3), -TotOut (1, 5) AED 


velx=(0 0]; vely=(7000 7000]; velz=(4500 1.S*Misealey (1) 44500) 
plot3(velx,vely,velz,'r.-') 
ncyx-[0 0]; псуу=[6500 6500]; ncvz-2[4500 7*abs(omegaout(1,1))-74500]; 
plotsS(ncix псу поља n gn и) 
nezx=(0 0]; nezy=(6000 6000); nezz=(4500 7*abs (omegaout(1,2))+4500]; 
pDliosjt4nczxneczv,ncez2 b 
TIN Emaxe-rounci tms и) 
plot3([MissileOut(].l1) TgtOut(l,1)]2[MySS ете (172% 
TgtoutüW»),... 
[-MissileOut(1,3) -TgtOut(1,5)],'k") 
end 


Rolda Ont 
grid on 


title(['Missile Engagement ',date],'FontSize',18) 
text (35000760007 6000) [num2stritmax, $22.21 7 ИЛ 
seconds'], 'FontSize',18) 
text (0; 7000,1.5*MissileV(1)+4500 1 
num2str(MissileV(1); %4.0f/)], FontSize’, 14) 
text(0,6500,7*abs(omegaout(1,1))+4500,[ ” ' 
rnum2 str (omegaout (1, 1) Z79 8045 '%2 Е) "Более ІШІ) 
text(0,6000,7*abs(omegaout(1,2))+4500, [ ' í 
num2stcr (omegaout (1, 2)79.8045,?%2.1f')];, PoncSize T3 


M(timestep) =getframe; 
end 


movie (M) 


%//end of file spielberg.m 
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funcELOn y=tgo(u) 


STGO Computes time to go from Range and Range Rate 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


qb f koh ke che ke coke eoe e e e e ke e ke cce ke ke ke ecce e heck kk ko kk eec ke ck ke ce e KR RK KK KK KK RK KK KR RK KK KK KKK К 
ЖЖЖ ЖЖЖЖ Ж 

$// File: tgo.m 

$// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

$// Date: 31 May 00 

$// Description: computes tgo 

$// Inputs: range, range rate 

#// Outputs: tgo 

$// Process: 

$// Assumptions: 

$// Warnings: 


$ / «ЖКХ 
Ххх хх xk 


%// Order of elements 


%// -Define globals 

$// -Define constants 

%// -Define elements of input vector 
$// -Functions 


I / É ХЕ ЖАЗСА АЖ Ж ЭАЖ ЭС ЖЖ ЖЭС ЖА ЖА ЖАТА СК КТК НЕ Ах Ааа 
kak kk kk KEK 


БЕ е хлле  сса 8°“ 
tac:  Gertine Constants  ****** 
и input Vector ^^ 


реак Panera li Zeevarlaples к кокк» 


% Kak KKK KK KK KE functions * * * * * * * * * * * К 
if (u(2)==0) 
x= 100; 
else 
y=abs (u(1)/u(2)); 
end 


$//end of file tgo.m 


[75 


funcErlonwstgtsetí(Range,Alt.Hdg) 
STGTSET Sets tgtinit variable for missile simulations 
% default tgt speed set at .83 Mach. Enter altitude 


% as a positive number 

% TGTSET (RANGE, ALT, HDG) 

% see also QUATERNION, BQUAT 

% Copyright 1999-2000 by Triple B Enterprises 


% / / Rok oe he oe ee e e he eee e e k Tk k Kk Kk k Kk kk kk k kk kk kk kk kk K kk 'k k e Se e he k * e c e e e k e e e kk kk 
ck oko ko ko ko ko kk 

$// File:  tgtset.m 

$// Name: LCDR Robert D. Broadston 

$// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 
$// Compiler: MatLab v5.3 

$// Date: 8 Jun 00 

$// Description:  initializes target state vector 

$// Inputs: see above 

$// Outputs: target state vector 

$// Process: 

$// Assumptions: 

$// Warnings: 


$5 f f NC e e ke ee e ke ke he ke ee e e ke he ke e e e e ke ke e e ehe ehe he echec he e e e he e ke e ke ke Se ke Sk ke e ke ke e e e ke e e ec ke ke Ж 
* * x xw xxx x 


%// Order of elements 


%// -Define globals 

%// -Define constants 

% / / -Define elements of input vector 
%// -Functions 


ужи лк КЛАВА КА А 
хжхххххххх 


55-5 о 


$ "define constants ****®=»* 
tgtmach=-. 83; % user sets Mach # 


qo ce черте три vector — ** 4224 


Ее а “ехе 
tgtspd-tgtmach*machvalt(Alt);$ machine computes speed 


% ххххххххкхххх functions ске хх к 
у= (Range; cos (Hdg*pi/180) *tgtspd; 0-sin(Hdg*pi/ 160) tese =o = ME 
% Note: negative altitude is for NED coords 


$//end of file tgtset.m 
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tunctryou vsthHebigstop(u) 
S&THEBIGSTOP consolidated simulation stop function 


% 
% see also 
% Copyright 1999-2000 by Triple B Enterprises 


if fe ааа а KKK AK KK KR RK KAR AR ARK А A K k kk k k XK SO NOE x x x 
ххххххххх 

%// File: thebigstop.m 

%// Name: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

$// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

$// Date: 17 Sep 00 

$// Description: stops simulation under variety of conditions 
%// Inputs: see below 

%// Outputs: -stop frag 

%// Process: 

%// Assumptions: 

%// Warnings: 


ШЕЕ Ж E E ое 
хжххххххх 


%// Order of elements 


% / / -Define globals 

%// -Define constants 

% / / -Define elements of input vector 
%// -Functions 


Е 
ххххххххх 

$ ****** define globals - ^- 

global stopflag 


b; tita A J Spin O Constanks ж” 


Беле Дере трип Vector --- === 
R= је. 

РЕ: 

МЕ (У), 

Vt=u (4); 

G=u (5) ; 

Ny=u (6) ; 

Nz=u (7); 

@тте=п (б): 


$ ****** Initialize variables  ****** 
stop=[]; 
y=0; 


Гү 


— — а "опна ва та + #54 ни 2 4 ӨМГҮЛІЗ 


9; * * * k kk k xk ck ko koX funcblomnis * k Kk Kk kx x xx xk xk x 


% check cases 
%1f (G>700) 


% ws] I 
% stop=’G stop’; 
$end 


if ((time>2.0)&(Vm<Vt) ) 
yz111; 
stop=’V stop’; 

end 


if ((time>2.0)&(Rdot>0) ) 
у=111; 
stops'Rdot stop'; 
end 


if (R«1e-6) 
у=111; 
stop=' R stop'; 
end 


if ((Nz>30*9.8045) | (Му>30%9.8045)) 
s= his. 
stop=’Cmd stop’; 

end 

БЕ StOprlag 
disp(stop) 

end 


$//end of file thebigstop.m 
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ЖУУК КККК Кккк хх ккк хх ххх ххх хх хх кк хх хх e hee ne e ecce e ккк КУ КККК хх хк 


ck cock ck ХХ x 


5// 
$// 
%// 
%// 
%// 
%// 
%// 
Б 
$// 
oye 
%// 
%// 
BA 


File:  thesis2plot.m 

Name:  LCDR Robert D. Broadston 

MSEE/EE Thesis 

Operating Environment: Windows NT 4.0 Service Pack 5 

Compiler: MatLab v5.3 

Date: 

Description: Plots a variety of data for missile tracking 
Simulations for use in thesis paper 

Topics: none 

Outputs: purdy pitchers 

Process: 

Assumptions: 

Warnings: 


БИ XC XOKCRUK OE KORR OR k k kk k k k K X K k N Kok kk oe oppo Ааа ат тт 


ххх ХХ ХХХ 


%// 
bu 
%// 
%// 
t 


Order of elements 
-Define globals 
-Define constants 
-Define elements of input vector 
-Functions 


Ф//ХХХХЕХХХХХХХХЕХХЕХХХЕХХХЕЖХЕКХХЕХХЕХЕХЕХХХХХХЕКТХЕКТХЕХХХЕХЕХХХХХЖЖЖХЖХХХЯЭ 


* * * w k xk ko kx ox 


и а Clo pels кај = 


mo ****** define Constants *****»* 


Б «ххх define Input vector Irita 


$ ****** Initialize variables ****** 
time=rem(now,1); 

= со Ее 24): 

пае= ови (кет( вале 24,1) 507. 

Емеса Киа сЕ ПЕ росе mne) 


% ЖЖЖЖ КЖК ЖЖ fare ls Om s KKKKKKKKK KEK 


$ engagement geometry 


figure(1) 

subpprlott2 5 1) 
plot(TgtOut(:,1),TgtOut(:,3),':',MissileOnt(:,]),M3ssembPeOUb 6 
axis equal 

outtextl=[’time: ’,num2str(ip),’ seconds’ ]; 

outtext2-['range: ',num2str(min(range)),' meters']; 


text (300,4000, ‘Intercept at:’) 
exe (300, 3500, outcexc. } 
text (300, 3000,outtCtext2) 


title(’Engagement Geometry’ ) 
xlabel('x (meters) ’) 
ylabel(’y (meters) ’) 
legend(’Target’, ’Missile’) 


1715 


en #4 OF 24 SE 323 FLL 


% missile to target distance 
range-sqrt((MissileOut(:,1)-TgtOut(:,1)).^2-*... 


(MissileOut(:,2)-TgtOut(:,3)).^2*(MissxleOutQN pu On E 


t=MissileOut(:,14); 
t_disc=0:FILTSAMP:max (L); 
index=find(range==min(range)); 
ip=t (index(1)); 


subplot 2,172) 

plot(t,MissileV) 

title('Missile to Target Range') 

xlabel(['time (seconds) ',date,timestr]) 
ylabel('missile velocity (m/s)"') 


$ missile accelerations 

gtorcessqrt(AccelOut(:w).^2-AÀccelOut( wm) TP 
"Ассе Оши: > )-72).7/9.80245- 

Ејацхке (2) 

suüubplot(2,1,1) 

plot Ee c EGO) 

title('Missile Accelerations') 


ylabel (‘Acceleration (g)’) 

axis([0 round(max(t)) 0 50]) 

6 сопршсе cose Lumetion J=20*e(tL) “2+integiuZ)7200 

u2= (omegaout(:,1).*2+omegaout(:,2).%2); 

їпредга!=О: 

ТӨК 11-221 ех 
integral=integral+(t(ii)-t(ii-1))*u2(ii-1); 


end 

J=20*min (range) *2+integral/1000; 

endfor 

outtxts['MQissile drvert:* ",num2str(integral)]; 


xlabel (outtxt) 
$ guidance command 
= Баје оде КА 1:2) 


plot (t omegaout(:, 1), c өтедесше 2). о) 
title(’Guidance law command output’) 

OUEExG=|"EOSt У: “ВЕ о time (seconds)']; 
xlabel([outtsxt;: ', date,timestr]) 


ylabel('n c (m/sec^2)") 
axis([0 round(max(t)) -300 3001) 
legend “<n e y у n ez’) 


$//end of file thesis2plot.m 
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2) 


DP еее жж тж а к жж ЖЖЖ Эх A XK 
Kk * * * ХХ * 


$// File: thesisinit.m 

%// Мапе: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 
%// Compiler: MatLab v5.3 

%// Date: 5 Sep 00 

%// Description: This script file initializes thesis work 
$// missile simulation 

$// Inputs: none 

$// Outputs: none 

%// Process: 

%// Assumptions: 


%// Warnings: 
S kk ck ok kk E OE EORR REO RE та ааа атта аа тука аа а тата аа та ek Te: 


ok ook kkk KK 


%// Order of elements 


%// -Define globals 

% / / -Define constants 

$// -Define elements of input vector 
$// -Functions 


6/ G KAKEK A S K X e e ok KC ke oie EX OK X X Z X X X ZX X X X X X X X K X X X X X ck C oe ok oe oko oe ck ok che che Ko ХАЖ жж 


*ck ko ko kock k k x 


clear 
бен цер түш а стт ккк 
global stopflag satflag XLAST FILTSAMP 


Ronee es) Шертпе сопесашсе и 

% physical constants 

omega_x=7.292115e-5; % earth rotation rate 

ӨМ Е-3.9860014е14; % G*mass of earth 

Ү Е-6.378164е6: % radius of earth 

r= J 205 2575 % ellipsoidal squash factor 
omega_E=[omega_x;0;0]; $ earth rotational velocity vector 


Ecce Celine inoue Vectors xs 


go oreka eitia] ze variables ии я 
$ missile physical parameters 
MissileData; 

% ахачлт:551]е; 


% missi e velocity vector 
vD |270: 0:01; 


% initial missile position vector 
c= (070, 5000) *- $ note altitude is negative in NED coord 
tr 


$ compute Euler angles 
вел = О“ ра 180; 
theta-0*p1/180; 

Ola Opa 19D. 


18] 


9; Kk Kk k k kx x x xxx x fs ll ei IG) rara *k * K Kk xx xx xxx x 


q O=quaternion(phi,theta,psi); 
q-0sqocvysart(g 0(1)^24920(2) 2*0 .0(3) ^2*Q9 ома не 


Bzquat2b(q 0); 


Р=0 раев 
БЕЛО SU 
во, 


omega_B=[P;Q;R]; 


% initial state vector 
init-[p;vB;omega B;q 01]; 


$ target initial state vector 
ЕЕЕ SDODS 

-250; 

OF 

250- 

-6000; 

o 


tmax=200; 
target turn-0; $ set target turn rate, default-0 
satflag-1; $ enable saturation of cmd accel 
ХЕҺЕӘШ- |25000; % initialize filter 

-250; 

Ü: 

OF 

250; 


FILTSAMP=.1; $ set filter sample interval 


$//end of file thesisinit.m 
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function y=vcpropnavpt (u) 
SPROPNAVPT Computes exact proportional navigation 


% with dragforce inputs for point mass simulation 
% see also EXACTPROPNAV2 
% Copyright 1999-2000 by Triple B Enterprises 


QS / f ok ck eoe ke oe oe coke oe X X XX KIKIY S XA К хк Кк КК ххх хх хх хк ххх ххх хх у ххх хк ЖЫКЫ КЕКЕК КЕКЕК X x x x 
хх Xx k k k x x 


$// File: VCPROPNAVPT.m 

%// Мапе: LCDR Robert D. Broadston 

%// MSEE/EE Thesis 

%// Operating Environment: Windows NT 4.0 Service Pack 5 

%// Compiler: MatLab v5.3 

%// Date: 24 May 2000 

%// Description: Proportional navigation guidance law for 6DOF 


%// flight model. Computes applied forces for use 
% / / by induced drag model. Required to eliminate 
% / / algebraic loops in the simulation 

%// Inputs: [seeker data, IMU data,timer] 

$// OUtDuts: [command accelerations,applied forces] 


$// Process: proportional navigation law 
$// Assumptions: none 
$// Warnings: none 


% / / * * * * ** k * he e ehe ee e ehe hee ee ee eee e e eee he e e ee e e hehe cec e e e ce e he e eoe ee e koe e e e e kc kc kk К Ж 
ххххххххх 


%// Order of elements 


% / / -Define globals 

% / / -Define constants 

%// -Define elements of input vector 
%// -Functions 


E EE KAAI А ИАТА РА 
ххххххххх 

БК еше оо *~*=*** 

global m satflag 


$ rr define Constants “****** 
Nprime=5; 
Nprimez=5; 


Б М define трос мес ктк» 
thetadot=u(1); 
phidotzu(2); 
lossu(3); 
phrtossuud 
ВЕ 
Vc=-u(5) ; 
heading=u(7); 
МЕ 0507; 
Vmdot=u (9) ; 
=, 
Ебрева=ш 5) 
psizu(12); 
змеи 
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% ****** initialize variables  ****** 


% k k Kk xx wx x xx x x tunctiors ххжххххххххх 


ny=Nprime*750* (thetadot) /cos(psi-los) -Vmdot*tan(psi-los); 
nz=Nprimez*750* (phidot) /cos(theta-philos) -Vmdot*tan(theta-los) -9.8045; 


% control force limiter 
if satflag 
if (abs(ny)»30*9.8045) 
Av—-cugmi may) У Оле d 
епа 
Е 
nz-sign(nz) 30:9. 8045; 
end 
end 


% compute ABC forces applied 
Ex 
Fy-ny*m; 


FZ=nzZ Ш. 


$ Output vector 
vsLmunz FEX Ey, FZ]; 


%//end of file PROPNAVPT.m 
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APPENDIX C. SIMULATION DATA 


This appendix contains the plots listed below for each guidance law. The 
engagement geometry is the same for each guidance law, initial range; 20 km, attack 


azimuth 45 degrees, co-altitude; 6,000 meters, 6 g target maneuver at 3 seconds fy. 


l. Plan view of engagement 

D Missile velocity profile 

3i Missile accelerations and target acceleration estimates for filtered laws 
4. Guidance law command accelerations 


The guidance laws are: 

1. PN with N=5 

2. VCPN with constant gain 
25 Bang-bang 

4. Differential games 

о: APN with A=5 

6. Noisy seeker, PN with N 25 


This appendix also contains plots from the full aerodynamic model running in an 


open loop with control surface deflections as the control input 
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VCPN WITH CONSTANT GAIN 


Engagement Geometry 
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BANG-BANG 


Engagement Geometry 
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D. 


DIFFERENTIAL GAMES 


Engagement Geometry 
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APN WITH A=5 


Engagement Geometry 
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Е. 


NOISY SEEKER, PN (N’=5) 
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С. 
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FULL AERODYNAMIC MODEL 
Missile Altitude vs. Time 
2.25 degree up elevator, 1000 m/s initial velocity 








6600 
Missile operating in open loop | 
control with steady state input \ 
on elevator deflection angle 
6400 | | 
|| ] 
6200 | 
6000 
5800 
5600 
time (sec) 
Angle of Attack Response 
2.25 degree up elevator, 1000 m/s initial velocity 
1 ЛЕР EE 
Missile operating in open loop | 
control with steady state input | 
_ on elevator deflection angle 
D 
Ф 
C» 
Ф 
= 
< 
O 
< 








4 


5 
time (sec) 


192 


acceleration (g) 
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